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SYMBOLS 


A j 

B 

C i 

C K 

CF i 

CM t 

D F n ,C } 

D U 

D ^j 

D F. r Pj 

°h A i 

°M i , Uj 

D Mj,Pj 

F‘(X) 
F l ( Y) 

Fc i 

FDUj 

FDPj 

FRj 


angle of attack 

linear acceleration component 
angle of sideslip 
direction cosine rate parameter 
Ath control surface 
/ th control force 
/ th control moment 

derivative of the nth aerodynamic force component with respect to the Ath control 
surface 

direction cosine 

derivative of the /th aerodynamic force component with respect to the /th component 
of linear velocity 

derivative of the /th aerodynamic force component with respect to the /th component 
of angular velocity 

derivative of the /th aerodynamic force component with respect to the /th component 
of linear acceleration 

derivative of the /th aerodynamic moment component with respect to the /th compo- 
nent of linear velocity 

derivative of the /th aerodynamic moment component with respect to the /th compo- 
nent of angular velocity 

force component in wind-tunnel axes 

force component in body axes 

/th component of gravitational force 

/th component of aerodynamic force due to linear velocity perturbations 
/th component of aerodynamic force due to angular velocity perturbations 
/th component of the inertia force 

iii 



/ th component of thrust 


FT i 


g 

H 


IMj 




L in 


MDU i 

MDP t 


[*i], 

[/? 2 ],[/? 3 J 


gravity vector 

angular momentum vector 

/th component of the inertia moment 

a component of the inertia tensor 

angle of elevation of the nth thrust vector 

/th component of the position vector of the point of application of the nth thrust 
vector 

/th component of the aerodynamic moment due to linear velocity perturbations 
/th component of the aerodynamic moment due to angular velocity perturbations 
/th component of the angular velocity vector 
rotation matrices 


SF i 



SMj 


X 

TD F t U, 

TD I,C 

TD F .p. 
r 1 


/th component of the static aerodynamic force in body axes 

nth component of static aerodynamic force in wind-tunnel axes 

/th component of static aerodynamic moment in body axes 

nth component of static aerodynamic moment in wind-tunnel axes 

transformed aerodynamic stability derivative of the / th aerodynamic force component 
with respect to the /th component of linear velocity 

/th component of the transformed control derivative 

transformed aerodynamic stability derivative of the / th aerodynamic force component 
with respect to the /th component of angular velocity 


TDp £ transformed aerodynamic stability derivative of the / th aerodynamic force component 
1 1 with respect to the /th component of linear acceleration 


TDjy y . transformed aerodynamic stability derivative of the /th aerodynamic moment compo- 

1 1 nent with respect to the /th component of linear velocity 


TD^j p transformed aerodynamic stability derivative of the /th aerodynamic moment compo- 
1 ! nent with respect to the /th component of angular velocity 


iv 



TMj / th component of the thrust moment 

Uj /th component of the linear velocity vector 

V velocity vector 

X' wind-tunnel coordinates 

Y' body coordinates 

Ay4 ; - linear acceleration increments 

A Cf( control increments 

A Pj angular velocity increments 

A Uj linear velocity increments 

cc angular velocity vector 

Subscripts 

E earth-fixed coordinate system 

F forces 

i components of aerodynamic forces and moments, gravity forces, inertia forces, thrust 

forces and moments 

M moments 

n force and moment components 

0 initial values 

Superscripts 

a identifying indices and indices of contravariance 

1 3 identifying indices and indices of contravariance 

1 identifying indices and indices of contravariance 

j identifying indices and indices of contravariance 
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COMPUTER FORMULATIONS OF AIRCRAFT MODELS FOR SIMULATION STUDIES 


James C. Howard 
Ames Research Center 


SUMMARY 


A mathematical model of a dynamical system describes the physical characteristics of the 
system and can be used to determine the response of the system to the forces encountered. In the 
case of aeronautical systems exposed to gravity, inertia, aerodynamic, and thrust forces, the mathe- 
matical model enables the engineer to determine the state of the vehicle, its spatial orientation, and 
the geographical location as a function of time. Sometimes the vehicle response is determined by 
direct calculation; at other times, simulators are used. The simulator moves in response to the com- 
puted solutions of the mathematical model equations. These solutions represent the response of the 
system to the forcing functions generated by the simulator pilot’s control inputs. The formulation 
of mathematical models of complex multidegrees-of-freedom dynamical systems is time consuming 
and subject to human error. In view of the complexity involved, it is desirable to mechanize as 
much of the formulation as possible. Recent developments in formula manipulation compilers have 
led to an extension of the area of application of digital computers beyond the purely numerical data 
processing stage. These developments, combined with the design of several symbol manipulation 
languages, enable computers to be used for symbolic mathematical computation. This technique 
provides for the symbolic manipulation of mathematical expressions; for example, the expression 
SIN(X) can be differentiated, resulting in the expression COS(X). Moreover, it can be used fre- 
quently to obtain symbolic solutions in problem areas that heretofore could only be approached 
numerically. A computer system and language that can be used to perform symbolic manipulations 
in an interactive inode have been used to formulate a mathematical model of an aeronautical 
system. The example demonstrates that once the procedure is established, the formulation and 
modification of models for simulation studies can be reduced to a series of routine computer 
operations. 


INTRODUCTION 


Excluding the control loops, the mathematical description of an aeronautical system requires 
at least 12 equations (ref. 1): 3 force equations; 3 moment equations; 3 Euler angle equations or 
9 direction cosine equations to determine the spatial orientation of the vehicle, and 3 equations to 
determine the geographical location of the vehicle in inertial space. In view of this complexity, it 
is desirable to mechanize as much of the formulation as possible. 

Research undertaken with the object of mechanizing the formulation of mathematical models 
of aeronautical systems, has directed attention to the use of digital computers for symbolic mathe- 
matical computation. To date, the majority of computer systems and languages has been developed 
to facilitate the processing of numerical data in one form or another. More recently, however, a 
variety of symbol or formula-manipulation languages has evolved. The choice of system and lan- 
guage to be used for a given purpose depends on accessibility, personal preference, the type and 



magnitude of the problems to be formulated, and the computer facilities available to the user. At 
the time of writing, the two most important contenders in the symbol manipulation field appeared 
to the author to be REDUCE and MACSYMA. REDUCE is a program designed for general algebraic 
computations of interest to mathematicians, physicists, and engineers. In addition to the usual alge- 
braic manipulations, it has the capability of performing calculations of special interest to high 
energy physicists. Originally, it began as a system for solving special problems that arise in high 
energy physics, where much tedious repetitive calculation is involved. However, it was quickly 
recognized that the computer processes being used were quite general, and could be used for a great 
variety of algebraic manipulations. Although REDUCE can operate in a batch processing mode, it 
is intended primarily for interactive calculations in a time shared environment. Hence, it is com- 
mand oriented, rather than program oriented, since the result of a given command may be required 
before proceeding to the next step. REDUCE is available on most IBM 360 or 370 series computers, 
the DEC PDP-10, and the CDC 6400, 6500, 6600, and 7600 machines. At the time of writing, the 
MACSYMA system was available only at MIT through the Advanced Research and Project Agency 
(ARPA) Network (ref. 2). It is a large computer programming system, which can be used to perform 
symbolic as well as numerical mathematical operations. It was developed by the Mathlab group of 
project MAC’s Automatic Programming Division specifically for interactive use. In addition to 
manipulating algebraic expressions, the MACSYMA system can differentiate, integrate, take limits, 
solve systems of linear or polynomial equations and factor polynomials, expand functions, plot 
curves, and manipulate matrices. Moreover, it is continuously evolving to meet the needs of users. 
In view of its flexibility and diverse capabilities, MACSYMA was the system chosen to formulate 
mathematical models of aeronautical systems for simulation purposes. Subsequent experience with 
the system confirmed the wisdom of this choice. Of special interest to users is the facility with 
which MACSYMA can derive special forms from a more general formulation. The same feature 
permits users to introduce new sets of system parameters or modify existing ones as the occasion 
demands. It will be seen that a simple programming statement can be used to introduce a new set of 
inertia tensors, static aerodynamic coefficients, control force derivatives, aerodynamic stability 
derivatives, or thrust coefficients to meet model or design changes. An important aspect of the 
formulation of mathematical models of aeronautical systems for simulation and other purposes is 
the specification of the system of forces and moments. In aeronautical applications, the thrust and 
inertia forces and moments, and the gravity force can be formulated without difficulty; but the 
aerodynamic forces and moments require more detailed consideration. These are represented by the 
static forces and moments, the control forces, and the perturbation forces that depend on the aero- 
dynamic stability derivatives. These forces and moments have to be transformed from wind or wind- 
tunnel stability axes to aircraft body axes before the formulation can proceed. Although these for- 
mulations and transformations are not complicated, they are complex and unwieldy and are likely 
to contain errors when performed manually. The interactive capability, versatility, and simplicity 
of the MACSYMA system make it attractive to programmers and nonprogrammers alike. To illus- 
trate these aspects of the system, a mathematical model of an aeronautical system has been formu- 
lated and subjected to a series of modifications. 
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ANALYSIS 


Transformation Laws 

A necessary preliminary to the formulation of mathematical models of aeronautical systems is 
the transformation of static aerodynamic forces and moments, control force derivatives, and the 
aerodynamic stability derivatives from wind or wind-tunnel stability axes to aircraft body axes. It 
will be seen that whereas the static forces and moments obey the same transformation law as the 
system coordinates, the aerodynamic stability derivatives transform like the components of a mixed 
tensor, having one index of covariance and one index of contravariance (ref. 3). Moreover, because 
of the equivalence of covariant and contravariant transformations in orthogonal Cartesian systems 
of coordinates, the transformations can be treated as doubly covariant or doubly contravariant, if 
this simplifies the formulation. The rule for transforming static force coefficients from the X frame 
of reference (the wind-tunnel axes system) to the Y frame (the body axes system) is obtained as 
follows. 

Since all vectors, including the position vector of a point, obey the same transformation law, 
it follows that the force and moment vectors obey the same transformation law as the system coor- 
dinates; that is, if Sy n denotes a static aerodynamic force in the X frame of reference, and SF t 
denotes the corresponding transformed force in the Y reference frame, then since the transforma- 
tion law for coordinates is 


it follows that 


yi _ d Y l 
dX n 


9 Y' 

SF = 5 

dX n n 


( 1 ) 


(2) 


where Y - Y(X) remains to be specified. 

Likewise, if /)/.- denotes the nth control force derivative with respect to the Kth control 
surface, as measured in the X reference frame, and TDj ( . denotes the corresponding transformed 
derivative in the Y frame, then 


3 Y 1 

TD , = -- Du r 

I,C r..,L is 

dx" n K 


(3) 


The aerodynamic stability derivatives measure the rates of change of aerodynamic forces and 
moments with respect to motion vector components. In keeping with the usual practice in aerody- 
namic formulations, motion vector components will refer specifically to components of the linear 
velocity vector, components of the angular velocity vector, and components of the corresponding 
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linear and angular acceleration vectors. The transformation law for these derivatives may be 
obtained as follows: 

Let F l ( Y) be a force or moment in the Y system of axes, and let Ul(Y) be a 
motion vector component in this system of axes. 

Similarly, let F a (X ) be a force or moment in the X system of axes, and let 
IjP(X) be a motion vector component in this system of axes. 

Then the stability derivatives with respect to motion components, as measured in the Y system of 
axes, are related to the corresponding derivatives in the X system of axes by the following equation : 


dF‘( Y) = dF i ( Y ) bF a (X ) d(JP(X) (4) 

9 Ul(Y) dF a (X) dUP(X) dU>{ Y) 

Force, moment, and motion vector components obey the same transformation law as the system 
coordinates, that is, since 


Y i = an' x<x 
dX a 


F i (Y) = — F a (X ) 

ax a 


and 


ui{Y) = 


d XL 

dX@ 


U$(X) 


it follows that 


UP{X) = U>{Y) 

a y> 


Substitution from equations (5) and (7) in equation (4) gives 


dF‘(Y) = /3 Y‘ dX$\ dF a (X) 
3 Ul{Y) \?)X a dYl / dUP(X) 


(5) 


( 6 ) 


(7) 


(8) 


4 



In these and subsequent equations the summation convention is assumed; that is, if in any 
term an index oecurs twice, the term is to be summed with respect to that index for all admissible 
values of the index. 

Equation (8) shows that the aerodynamic stability derivatives transform like the components 
of a mixed tensor, having one index of covariance and one index of contravariance. Being a tensor 
of rank 2, equation (8) represents 9 equations, with each equation having, in general, 9 terms. 

A disadvantage of the transformation law as formulated in equation (8) is the requirement that 
the transformation equation X = X(Y) must be used. Fortunately, it is possible to avoid the use of 
the inverse transformation X = X(Y) if the coordinate transformations are orthogonal Cartesian. It 
can be shown that for orthogonal Cartesian transformations 

3 Y‘ a X> 

— 7 (9) 

a xJ a y 1 

Substitution of this relationship in equation (8) yields 

dF l (Y) = idY^ 3_W\ ZF a (X± (10) 

a u’{Y) \a^ a a X&) duP(X) 

The form of equation (10) shows that if the transformations are orthogonal Cartesian, the 
aerodynamic stability derivatives can be treated as doubly contravariant tensors. This form has the 
advantage that the inverse transformation X = 2f(T) is no longer required. 

Although the notation used in equations (8) and (10) reveals the tensor character of the 
transformation, superscripts will not be used in subsequent work. They will be replaced by sub- 
scripts, which are more convenient for programming purposes. For the same reason, the Greek sym- 
bols a and (3 will be replaced by the letters M and N. The following definitions are required: 

8 F a (X) _ dF ot( X) = _ D (] 1} 

a uP(X) Wp(X) dU N (X) f M’ u n 


and 


a f‘(Y) = dF j( Y ) = m 

a ui{Y) Wj(Y) F i’ u i 


( 12 ) 


where Uj\j denotes the derivative of the A/th component of aerodynamic force, with respect to 

the A/th component of the motion vector; and TDp. \j. denotes the transformed derivative of the 
;th component of aerodynamic force, with respect ^to^ the y'th component of the motion vector. 
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Aeronautical Reference Systems 


There are many coordinate systems in use in aeronautical research. Aerodynamic data obtained 
from wind-tunnel experiments may be referred to wind axes or to wind-tunnel stability axes. When 
the wind axes are used, the X 2 axis is aligned with the relative wind at all times. Most wind-tunnel 
data are referred to the wind-tunnel stability axes system. For this system, the X\ axis is in the 
same horizontal plane as the relative wind at all times. In addition to the wind axes and the wind- 
tunnel stability axes, there are other systems of axes fixed in the body and moving with the body. 
These are referred to as body axes. In aerospace applications, a body axis system has the Y 3 axis 
fixed along the longitudinal centerline of the body, the Y 2 axis normal to the plane of symmetry, 
and the Y 3 axis in the plane of symmetry. The equations of motion of aerospace vehicles are formu- 
lated with respect to body axes. The main advantage of these axes in motion calculations is that 
vehicle moments and products of inertia about the axes are constants. When the body axes are 
chosen so that the products of inertia vanish, they are known as principal axes. A system of axes, 
which is frequently used to study the stability of aircraft in the presence of disturbing forces that 
produce small perturbations, is the flight stability axes. This is an orthogonal system fixed to the 
vehicle, the Y j axis of which is aligned with the relative wind vector when the vehicle is in a 
steady-state condition, but then rotates with the vehicle after a disturbance as the vehicle changes 
angle of attack and sideslip. Some of these axes are shown in figure 1 (ref. 4). 




Figure 1.- Systems of reference axes, including body, principal, wind, flight stability, and wind-tunnel stability. 


Transformation Equations 

The elements of the matrices defining a transformation from wind or wind-tunnel stability 
axes to body axes are functions of the angle of attack (A) and the angle of sideslip (B). Moreover, 
coordinates in wind-tunnel axes are denoted by a column vector of coordinates X/, and the body 
axes coordinates by a column vector Y/. To bring a reference frame from the wind axes into coin- 
cidence with the body axes involves a negative rotation {B) about the Y 3 axis, followed by a posi- 
tive rotation ( A ) about the Y 2 axis. These matrices may be entered and multiplied when communi- 
cation has been established and the system prints (C'l). When this occurs, the user types 
ENTERMATRIX(m.n) which allows one to enter a matrix, element by element, with MACSYMA 
requesting values for each of the (m,n) entries as follows: 
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(Cl) ENTERMATRIX(3,3); 


ROW 

1 

COLUMN 

1 

COS(A) ; 

ROW 

1 

COLUMN 

2 

0; 

ROW 

1 

COLUMN 

3 

-SIN(A) ; 

ROW 

2 

COLUMN 

1 

0; 

ROW 

2 

COLUMN 

2 

1; 

ROW 

2 

COLUMN 

3 

0; 

ROW 

3 

COLUMN 

1 

SIN(A); 

ROW 

3 

COLUMN 

2 

0; 

ROW 

3 

COLUMN 

3 

COS(A) ; 

MATRIX-ENTERED 



(Dl) 

(C2) ENTERMATRIX( 3,3) ; 


COS(A) 0 
0 1 
SIN(A) 0 


ROW 

1 

COLUMN 

1 

COS(B); 

ROW 

1 

COLUMN 

2 

-SIN(B) ; 

ROW 

1 

COLUMN 

3 

0; 

ROW 

2 

COLUMN 

1 

SIN(B); 

ROW 

2 

COLUMN 

2 

COS(B); 

ROW 

2 

COLUMN 

3 

0; 

ROW 

3 

COLUMN 

1 

0; 

ROW 

3 

COLUMN 

2 

0; 

ROW 

3 

COLUMN 

3 

1; 


- SIN(A) ] 

] 

0 ] 

] 

COS (A) ] 



MATRIX-ENTERED 


(D2) 


(C3) ENTERMATRIX(3,1 ) ; 

ROW I COLUMN 1 X[1 ]; 
ROW 2 COLUMN 1 X[2] ; 
ROW 3 COLUMN 1 X[3]; 
MATRIX-ENTERED 


(03) 


[ COS ( B ) 

- SIN(B) 

0 ] 

C 


] 

[ SIN(B) 

COS(B) 

0 ] 

c 


1 

[ 0 

0 

1 ] 


[ X 
[ 

[ 

[ X 
[ 

[ 

[ X 


] 

] 

] 

] 

] 

] 

] 

] 


(C4) (D1).(D2).(D3); 


[ COS(A) (X, COS (B ) - X 9 SIN(B)) - X 9 SIN(A) ] 
[12 3 ] 

[ ] 

(D4) [ X, SIN(B) + X 9 COS(B) | 

[ 1 2 ] 

[ ] 

[ SIN(A) (X-, COS(B) - X 9 SIN(B) ) + X 9 COS(A) ] 

[12 3 ] 

( C 5 ) FOR 1:1 THRU 3 DO ROW [ 1 ] : F I RST ( ROW ( ( D4 ) , I ) ) $ 

(C 6 ) FOR 1:1 THRU 3 DO (Y[I]:R0W[I][1],DISPLAY(Y[I])); 


Y 1 = COS (A) (X ] COS(B) - X 2 SIN(B) ) - X 3 SIN(A) 
(D 6 ) Y 2 = X ] SIN(B) + X 2 COS(B) 

Y 3 = SIN(A) (X ] COS(B) - X 2 SIN(B) ) + X 3 COS(A) 



In order to more fully appreciate the results obtained so far, the reader should note that 
MACSYMA requests the i th row and the ;'th column of the matrix being entered by typing 
ROWICOLUMNJ. The user merely provides the corresponding element. When all m X n elements 
have been entered, the system types MATRIX-ENTERED, formulates the matrix and assigns an 
identifying number (DI). When the user types the command (C4), that is, (D1).(D2).(D3), the three 
matrices are multiplied in the order requested and the product matrix is displayed in (D4). 

The two programming steps shown in (C5) and (C6) lead to the functional form (D6), which 
represents the required transformation from wind axes Ay to body axes Y j. 


FORCES 


Transformation of Static Forces 

The static aerodynamic forces transform like the components of a contravariant vector; that is, 
if Sp denotes a static aerodynamic force in the X frame of reference, and SFf denotes the corre- 
sponding transformed force in the Y reference frame, then from equation (2) 


SF- = u '- Sr? (2) 

bX n n 

where Y = Y(X) is obtained from the displayed output (D6). 

Given the transformation equations (D6), the transformed aerodynamic static forces are 
obtained by expanding equation (2). Three programming steps are sufficient to formulate the 
required values. The simple program and the displayed results are 

( C7 ) SF[I]:=0$ 

(C8) FOR 1:1 THRU 3 DO FOR M : 1 THRU 3 DO 
SF[I]:SF[I]+DIFF(Y[I],X[M])*S[F[M]]$ 

(C9) FOR 1:1 THRU 3 DO DISPLAY (SFC I ] ) ; 


SF, = - S F COS(A) SIN(B) + S F COS(A) COS(B) - S F SIN(A) 

i i- 2 h i 

SF- = S F SIN(B) + S F COS(B) 

6 1 y 2 

SF- = - S F SIN(A) SIN(B) + S F SIN(A) COS(B) + S F COS(A) 
6 2 h l h 3 

(D9) DONE 
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Transformation of Control Force Derivatives 


The control force derivatives obey the same transformation law as the static forces; that is, if 
Dy n Cj( denotes the nth control force derivative with respect to the Aith control surface as mea- 
sured in the X reference frame, and TDj (j denotes the corresponding transformed derivative in the 
Y frame, then using equation (3) 


7)Y l 

TDj „ = — — Dp r 

’ dx n Fn,Ck 


(3) 


where Y - Y(X) is again obtained from the displayed output (D6). 


As in the preceding section, the transformed control derivatives are obtained by expanding the 
transformation law for derivatives given the transformation equations (D6). The transformed deriva- 
tives are obtained by executing the following simple program, which has exactly the same form as 
the program used to transform the static forces. These are: 

(CIO) TD [ I ,C] : =0$ 

(Cll) FOR 1:1 THRU 3 DO FOR M:1 THRU 3 DO 
TD [ I , C ] : TD [ I ,C]+DIFF( Y[ I] ,X[M] )*D[F[M] ,C[K]]$ 

(Cl 2) FOR 1:1 THRU 3 DO DISPLAY(TD[I ,C]) ; 


TD 


r = -D F r C0S(A)SIN(B) + D F r C0S(A)C0S(B) - D F r SIN(A) 
i,l r 2 ,t K V K h 3 ,L K 


TD 3,C 

(D12) 


TD ? r = Dp r SIN(B) + Dp r COS(B) 

' 1 5 ^ K ' 2 9 ^ K 

-Dp r SIN(A)SIN(B) + Dp r SIN(A)C0S(B) + D r r COS(A) 
r 2 ,c K r ]? L K r 3 ,L K 

DONE 


The corresponding control forces are obtained by multiplying the control derivatives by the 
appropriate control increments A Cj^. The following two programming steps are sufficient to ensure 
evaluation of the required forces. These are denoted by CF [ in the displayed output. 

(Cl 3) FOR 1:1 THRU 3 DO CF[I ] :TD[ I ,C]*DEL(C[K] )$ 

(Cl 4) FOR 1:1 THRU 3 DO D I S PLA Y ( C F [ I ]) ; 


CF-, = (-Dp r COS (A)SIN(B ) 
r 2 ,L K 


+ Dp r C0S(A)C0S(B) - Dp r SIN(A)) DEL(C 1/ ) 
1 ,l K ' 3 st K K 
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CF 2 = (Dp c SIN(B) + Dp c COS(B))DEL(C k ) 
1 5 K 2 ’ K 


CF 3 °F, 


r SIN(A)SIN(B ) + D P r SIN(A)COS(B) + D F r COS(A)) DEL(C k ,) 

,l k i- 1 ,t |< i- 3 ,l k 


(D14) DONE 


Forces Produced by Linear Velocity Perturbations 

The next step in the formulation involves the determination of the aerodynamic forces pro- 
duced when an aircraft is subjected to linear velocity perturbations A U;. Before these forces can be 
determined, the aerodynamic stability derivatives, with respect to linear velocity components, must 
be transformed from wind or wind-tunnel stability axes to aircraft axes. For a detailed discussion of 
the transformation of these derivatives, the reader is referred to equations (4) through (10). In this 
application, the aerodynamic stability derivatives of the / th force with respect to the /th velocity 
components are denoted by Dp. p-. The corresponding transformed derivatives are denoted by 
TDp. {/.. The program required for this application assumes the form 

h J 

(Cl 5 ) TDU[I , J ] : =0$ 

(Cl 6) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
FOR M: 1 THRU 3 DO FOR N : 1 THRU 3 DO 

TDU[I,J]:TDU[I,J]+DIFF(Y[I],X[M])*DIFF(Y[J],X[N])*D[F[M],U[N]]$ 

It only remains to multiply the transformed derivatives by the appropriate velocity increments 
to obtain the required forces, which are denoted by FDU[. The next three programming steps 
instruct MACSYMA to evaluate and display the forces produced by linear velocity perturbations. 
These are 


( C 1 7 ) FDU[I ] : =0$ 

( Cl 8 ) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
FDU[I]:FDU[I]+TDU[I,J]*DEL(U[J])$ 

( C 1 9 ) FOR 1:1 THRU 3 DO DISPLAY( FDU[I ] ) ; 

FDU-J = (Dp^ C0S(A)SIN(A)SIN 2 (B) 

- D F C0S(A)SIN(A)C0S(B)SIN(B) 

r 2’ u l 

- Dp .. C0S(A)SIN(A)C0S(B)SIN(B) + D F „ S I N 2 ( A )S I N ( B ) 

1 ’2 h 3’ u 2 


1 1 



- D c . . COS 2 (A)SIN(B) + Dp „ COS(A)SIN(A)COS 2 (B) 

F 2’ U 3 1’ 1 

-Dr- 1 1 SIN 2 {A)COS(B) + Dp ,, C0S 2 (A)C0S(B) 

F 3’ u i *‘r u 3 

- D c COS(A)SIN(A) )DEL(Uo) 

F 3 ,U 3 


(-Dp „ COS(A)SIN 2 (B) - Dp „ COS(A)COS(B)SIN(B) 
‘ 2 * U 1 r 2» u 2 


D 


r l 


M COS(A)COS(B)SIN(B) - D r „ SIN(A)SIN(B) 
U 1 h 3’ U l 

D c .. COS(A)COS 2 (B) - Dp „ SIN(A)COS(B))DEL(U 
h l» u 2 3 ,u 2 

D r 


,) 


D c „ COS 2 (A)SIN 2 (B) - Dp „ COS 2 (A)COS(B)SIN(B) 

' 2 ’ U 2 r 2’ u l 

)p ,, C0S 2 (A)C0S(B)SIN(B) + Dp .. COS(A)SIN(A)SIN(B) 
F 1’ U 2 3 ,U 2 

) c . | COS(A)SIN(A)SIN(B) + Dp ,, COS 2 (A)COS 2 (B ) 

F 2’ U 3 r 1 

) c ,, COS(A)SIN(A)COS(B) - D F „ COS(A)SIN(A)COS(B ) 
F 3’ U 1 h l’ U 3 

)p ,, SIN 2 (A) ) DEL ( U ) 
r 3’ 3 

= (-Dp „ SIN(A)SIN 2 (B) - Dp ,, SIN(A)COS(B)SIN(B) 
f 1’ U 2 ' 2 ,u 2 


, SIN(A)COS(B )SIN(B ) + D F ,, COS(A)SIN(B) 
J 1 h l’ U 3 

. SIN(A)COS 2 (B) + D 
J i 


COS(A)COS(B))DEL(Uj 
2’ U 3 13 


F 7 ,U 


SIN 2 (B) + Dp . ■ COS(B ) S I N ( B ) + D F ,, COS(B)SIN(B) 
1 F 2’ U 1 1 ’ 2 

C0S 2 (B))DEL(U 2 ) + (-d f y COS(A)SIN 2 (B) 


COS( A)COS(B )SIN(B) + D 


COS(A)COS(B)SIN(B) 



F 1 ,U 3 


F 2’ U 3 


FDU 3 = (D 


F 1 ,U 2 


F 2 ,U 3 


F 3’ U 1 


F 3 ,U 3 


f 2 ,u 2 


F 3’ U 1 


F 3’ U 2 


F 2’ U 1 


F 1 ,U 2 


F 3 ,U 2 


F 1 ,U 3 


F 3’ U 3 


SIN(A)SIN(B) + D c „ COS(A)COS 2 (B) 
h 2’ U l 

SIN(A)COS(B))DEL(U 1 ) 


,, SIN 2 (A)SIN 2 (B) - Dp 1 1 SIN 2 (A)COS(B)SIN(B) 

SIN 2 (A)COS(B)SIN(B) - Dp .. COS(A)SIN(A)SIN(B ) 

1-3, u 2 

COS(A)SIN(A)SIN(B) + D p .. SIN 2 (A)COS 2 (B) 

*"1 * U 1 

COS(A)SIN(A)COS(B) + D F .. COS(A)SIN(A)COS(B) 

,u 3 

C0S 2 (A))DEL(U 3 ) + (-Dp u SIN(A)SIN 2 (B) 

SIN(A)COS(B)SIN(B) + D SIN(A)COS(B)SIN(B) 

h l ,U 1 

COS(A)SIN(B) + Dp n SIN(A)COS 2 (B) 

h ,U 2 

C0S(A)C0S(B))DEL(U 2 ) 


COS(A)SIN(A)SIN 2 (B) 

2 

COS(A)SIN(A)COS(B)SIN(B) 


COS(A)SIN(A)COS(B)SIN(B) + D F SIN 2 (A)SIN(B) 

' 2 ,u 3 

COS 2 (A)SIN(B) + Dp 1 1 COS(A)SIN(A)COS 2 (B) 

,U 1 

SIN 2 (A)COS(B) + Dp C0S 2 (A)C0S(B) 
h 3’ U l 

COS(A)SIN(A) ) DEL(U 1 ) 


13 



Forces Produced by Angular Velocity Perturbations 


The program used in the preceding section can, with suitable notational changes, be used to 
formulate the forces produced by angular velocity perturbations. However, whereas in the preceding 
application the required forces were obtained by multiplying the transformed aerodynamic stability 
derivatives by linear velocity increments, in the present case the transformed derivatives must be 
multiplied by angular velocity increments. In view of these similarities, the following program and 
displayed forces will be presented without further comment, except to point out that the aerody- 
namic stability derivatives of the / th force with respect to the y'th angular velocity component are 
denoted by Dp it Pj. The corresponding transformed derivatives are denoted by TDfiPh and the 
resulting forces by FDP[. 

(C20) TDP[I,J]:=0$ 

(C21) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
FOR M: 1 THRU 3 DO FOR N:1 THRU 3 DO 

TDP[I,J]:TDP[I,J]+DIFF(Y[I],X[M])*DIFF(Y[J],X[N])*D[F[M],P[N]]$ 

(C22) FDP [ I ]:=0$ 

(C23) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
FDP[ I ] : FDP [I ]+TDP[I ,J]*DEL(P[J])$ 

(C24) FOR 1:1 THRU 3 DO DISPLAY(FDP[I])'; 


FDP-| = (D p p COS(A)SIN(A)SIN 2 (B) 

- D f D COS (A)SIN(A) COS (B)SIN(B) 

K 1 

- Dp p COS (A)SIN(A) COS (B)SIN(B) + D P p S I N 2 ( A ) S I N ( B ) 

r 1 ’ k 2 h 3 ’ 2 

- Dp p C0S 2 (A)SIN(B) + Dp p C0S(A)SIN(A)C0S 2 (B) 

2’ r 3 r l ,r l 

- Dp SIN 2 (A)C0S(B) + Dp d COS 2 (A)COS(B) 

r 3’ i ^l ,l 3 

- Dp D COS (A) SIN (A) )DEL(P 0 ) 

h 3 ,P 3 

+ (-Dp p C0S(A)SIN 2 (B) - Dp p C0S(A)C0S(B)SIN(B) 
r 2 5 ' 1 * 2 * ' 2 

+ Dp p C0S(A)C0S(B )SIN(B ) - Dp p SIN(A)SIN(B) 
t ‘l ,p l 3 ’ K 1 


+ Dp D C0S(A)C0S 2 (B) - Dp p SIN(A)C0S(B))DEL(P 9 ) 
r 1 ’ k 2 3’ 2 



I 


+ (Dp p cos 2 (a)sin 2 (b) - d f p cos 2 (a)cos(b)sin(b) 

* 2 9 ' 2 r 2’ 

- Dp p C0S 2 (A)C0S(B)SIN(B) + Dp p COS(A)SIN(A)SIN(B) 

i ’ ' 2 . 3’ 2 

+ Dp p COS(A)SIN(A)SIN(B) + D F p COS 2 (A)COS 2 (B) 
h 2 ,l 3 i 

- Dp p COS(A)SIN(A)COS(B ) - D F D COS(A)SIN(A)COS(B) 

h 3’ P l i ‘r l 3 

+ D p p sin 2 (a))del(p 1 ) 


FDP 9 = (-Dp p SIN(A)SIN 2 (B) - D F p SIN(A)COS(B)SIN(E 
£ * 1 9 • 2 v 2' V 2 

+ Dp p SIN(A)COS(B)SIN(B) + D F p COS(A)SIN(B) 

M’ 1 ] l’ K 3 

+ d f p SIN(A)COS 2 (B) + d f p cos(a)cos(b))del(p 3 ) 

+ Dp D SIN 2 (B) + Dp COS(B)SIN(B) + D F p C0S(B)5 
h 1 ’ K 1 h 2’ h l h l ,f 2 

r. r\ 


+ Dp p COS 2 (B ) )DEL(P ? ) + (-Dp p C0S(A)SIN 2 (B ) 

' 2*' 2 r-| 5^2 

- Dp p COS(A)COS(B)SIN(B) + D F p COS(A)COS(B)SIN(B) 

r 2 ’ r 2 r l* 1 

- Dp D SIN(A)SIN(B) + Dp p COS(A)COS 2 (B) 

^1’ k 3 ^2 ,k 1 

- Dp p SIN(A)COS(B))DEL(P,) 

r 2’ 3 

FDPo = (Dp p SIN 2 (A)SIN 2 (B) - D F p SIN 2 (A)C0S(B )SIN(B ) 
j ' 2 5 ' 2 ' 2 9 ' 1 

- Dp p SIN 2 (A)COS(B)SIN(B) - Dp p COS(A)SIN(A)SIN(B ) 

v V v 2 3’ 2 

- Dp COS(A)SIN(A)SIN(B) + D F p SIN 2 (A)COS 2 (B ) 

' 2’ K 3 h l ,h l 

+ Dp p C0S(A)SIN(A)C0S(B ) + D F p C0S(A)SIN(A)C0S(B) 
h 3 ,K i 3 


Dp p C0S(A )SIN( A)SIN(B ) 
3 ,h 2 

ip D SIN 2 (A)COS 2 (B) 
h l ,K 1 
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+ D. D C0S 2 (A))DEL(PJ + (-Dp p SIN(A)SIN 2 (B) 

h 3 ,h 3 2’ 1 

- Dp D SIN(A)COS(B)SIN(B) + D F p SIN(A)COS(B)SIN(B) 

y 2' v 1 I s 1 

+ Dp p COS(A)SIN(B) + Dp p SIN( A)C0S 2 (B ) 

h 3 ,P l l’^2 

+ Dp p COS(A)COS(B) )DEL(P ? ) 

*"3*2 


+ 


+ 


+ 


+ 


(D 


f 2 ,p 


F P 

h 2’ K l 


F 1’ P 2 


T 2’ P 3 


F P 


T P 

' 3 ’ K 1 


COS(A)SIN(A)SIN 2 (B) 

2 

COS(A)SIN(A)COS(B)SIN(B) 


COS(A)SIN(A)COS(B )SIN(B ) 


SIN 2 (A)SIN(B ) - Dp p C0S 2 (A)SIN(B) 

3’ 2 

C0S(A)SIN(A)C0S 2 (B) - Dp p SIN 2 (A)C0S(B) 

^ ,h 3 

C0S 2 (A)C0S(B) - Dp p COS(A)SIN(A) ) DEL ( P ) 

*"3 ’ 3 1 


Forces Produced by Linear Acceleration Perturbations 

The procedure used in the preceding two sections may, with equal facility, be used to formu- 
late the aerodynamic forces produced by linear acceleration perturbations. However, in this case the 
required forces are obtained by multiplying the transformed aerodynamic stability derivatives, with 
respect to acceleration components, by linear acceleration increments. The aerodynamic stability 
derivatives of the i th force component Ff with respect to the /th linear acceleration component Aj 
are denoted by Dpj Aj, and the transformed derivatives by TDp^Aj. The corresponding force com- 
ponents in body axes are denoted by FDAj. 

Due to the fact that lift responds in a transient manner when, for example, the angle of attack 
A or the linear velocity component U 3 is suddenly changed, the acceleration derivatives are very 
different from the velocity derivatives, which can be determined on the basis of steady-state aero- 
dynamics. This is a consequence of the fact that the pressure distribution on a wing or tail surface 
does not adjust itself instantaneously to its equilibrium value when the angle of attack or the veloc- 
ity components are suddenly changed. Hence, in order to get a sufficiently accurate description of 
these derivatives during the indicial response phase, it may be necessary to use function generation 
or look-up tables (ref. 5). 
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When the program of the preceding section has been modified to incorporate the necessary 
notational changes, it assumes the following form: 


(C25) TDA[I,J]:=0$ 

(C26) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
FOR M: 1 THRU 3 DO FOR N:1 THRU 3 DO 

TDA[I,J]:TDA[I,J]+DIFF(Y[I],X[M])*DIFF(Y[J],X[N])*D[F[M],A[N]]$ 
(C27 ) FDA[I]:=0$ 

(C28) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
FDA[I ] : FDA[I]+TDA[I , J ]*DEL(A[J] )$ 

(C29) FOR 1:1 THRU 3 DO D ISPLAY ( FDA[ I ] ) $ 


Execution of this program yields the aerodynamic forces produced by linear acceleration 
:urbations. These are 

O 


per- 


FDA ] = (D F A COS(A)SIN(A)SIN 2 (B) 

. COS(A)SIN(A)COS(B)SIN(B) 
? »"1 


- D 


. C0S(A)SIN(A )C0S (B )SIN(B ) + D p . SIN 2 (A )SIN (B ) 
a 2 3 ,H 2 

D f . C0S 2 (A1SIN(B) + Dp . C0S(A)SIN(A)C0S 2 (B) 
h 2 ,A 3 1’ 1 


Dr 


. SIN 2 (A)C0S(B) + Dp . C0S 2 (A)C0S(B) 

b ,m i r 3 


- Dp . COS(A)SIN(A)DEL(A~) 
h 3 ,A 3 

+ (-Dp . C0S(A)SIN 2 (B) - Dp . C0S(A)C0S(B)SIN(B) 
r 2 » “ i 


+ D c 


'1 


. C0S(A)C0S(B)SIN(B ) - D F . SIN(A)SIN(B) 
’ A 1 h 3’ A l 


+ d f A C0S(A)C0S 2 (B) - Dp A ^ sin(a)cos(b))del(a 2 ) 

+ (Dp . C0S 2 (A)SIN 2 (B) - Dp . C0S 2 (A)C0S(B)SIN(B) 
r 2 » M 2 * Z 5 

- D c . C0S 2 (A)C0S(B)SIN(B) + Dp . COS(A)SIN(A)SIN(B ) 

r 1 5 r^2 ** 0 » 
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) F . COS(A)SIN(A)SIN(B) + D F . COS 2 (A)COS 2 (B ) 
h 2 ,A 3 1’ 1 

V . COS(A)SIN(A)COS(B) - D F . COS(A)SIN(A)COS(B) 
h 3’ A l 1’ 3 


D p A SIN^(A) )DEL(A-, ) 


^2 = ( ' D F r A 2 SIN(A)SIN 2 (B) - D p ^ 


SIN(A)COS(B ) S I N ( B ) 


+ D c . SIN(A)COS(B)SIN(B) + D F . COS(A)SIN(B) 

F 1’ A 1 ' 1 ,A 3 

+ D c . SIN(A)COS 2 (B) + D f . COS(A)COS(B) )DEL(Aq) 
h 2’ A l r 2 ,a 3 

+ (D c . SIN 2 (B) + D f . COS(B)SIN(B ) + D F . COS(B)SIN(B) 

r i , A-| V 2 9 \ ‘ 1 ’ M 2 

+ D f A C0S 2 (B))DEL(A 2 ) + (" D Fl ,A 2 COS(A)SIN 2 (B) 

- D c . COS(A)COS(B)SIN(B) + D F . COS(A)COS(B)SIN(B) 

r 2 5^2 ' 1 

- D f . SIN(A)SIN(B) + D f . COS(A)COS 2 (B) 

h r A 3 2’ 1 

- D f . SIN(A)COS(B ) DEL ( A ) 

F 2 ,A 3 


FDA., = (Dr- A SIN 2 (A)SIN 2 (B) - D F . SIN 2 (A)COS(B)SIN(B) 

o r 2 5 ^2 r ^ , h -| 


- D c . SIN 2 (A)C0S(B)SIN(B) - D f . COS(A)SIN(A)SIN(B) 

h l’ A 2 3’ 2 

- D f . COS(A)SIN(A)SIN(B ) + D F . SIN 2 (A)C0S 2 (B) 

f 2’ A 3 1’ 1 

+ D f . C0S(A)SIN(A)C0S(B ) + D F . C0S(A)SIN(A)C0S(B) 
F 3 ,A 1 1’ 3 

+ D f . C0S 2 (A) )DEL(Aq) + (-D f . S I N ( A ) S I N 2 ( B ) 

r 3’ 3 2 ,a 1 

- D F . SIN(A)C0S(B )SIN(B) + D F S I N ( A ) COS ( B ) S I N ( B ) 

r 2 5 ^2 ' 1 5 
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+ 


+ 


+ 


+ 


+ 


+ 


F 3 ,A 1 


F A 
r 3’ 2 


(D 




F A 
h 2 ,A 1 


J F,,A 


1 ’2 


F A 

h 2’ 3 


F A 
h l ,A 1 


F A 

h 3 ,A 1 


COS(A)SIN(B ) + Dp . SIN(A)C0S 2 (B) 

h l ,A 2 

cos(a)cos(b))del(a 2 ) 

COS(A)SIN(A)SIN 2 (B) 

cos(a)sin(a)cos(b)sin(b) 

COS(A)SIN(A)COS(B)SIN(B ) 

SIN 2 (A)SIN(B) - Dp . COS 2 (A)SIN(B) 

' 3 ,A 2 

COS(A)SIN(A)COS 2 (B) - Dp . SIN 2 (A)C0S(B) 

h l ,A 3 

C0S 2 (A)C0S(B) - Dp A COS(A)SIN(A))DEL(A 1 ) 


The components of the resultant aerodynamic force are 

(C30) FOR 1:1 THRU 3 DO FA[I ] : FDU[I ]+FDP[ I ]+FDA[I]+CF[I ]+SF[I]$ 


Gravity Forces 

The gravitational force vector acting on an aircraft has the value Mg, where M is the mass of 
the aircraft and g is the gravitational acceleration vector. The magnitude of g is assumed constant, 
which is tantamount to the assumption of a flat earth. The gravity vector is specified in an earth- 
fixed reference frame; and it is required to find the components of this vector in aircraft body axes. 
In accordance with aeronautical convention, a transformation from earth-fixed axes to aircraft body 
axes involves a rotation R 3 about the Y 3 body axis, followed by a rotation R 2 about the Y 2 body 
axis, and a rotation R i about the F, body axis. Hence, if it is assumed that the body axes and the 
earth-fixed axes are initially coincident, the components of the gravitational force FGj in body axes 
are given by an equation of the form 

[ FG ] = [R } |[* 2 ][/? 3 ][Afe] 

where (FG'J is a column vector of body axes components, [/?,], [/? 2 ] , and [R 3 J are rotation 
matrices, and [Mg j is a column vector of earth-fixed axes components. These matrix operations can 
be performed by MACSYMA to yield the required force components in body axes as follows; 
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( C 31 ) ENTERMATRIX(3,3) ; 


ROW 1 COLUMN 1 
ROW 1 COLUMN 2 
ROW 1 COLUMN 3 
ROW 2 COLUMN 1 
ROW 2 COLUMN 2 
ROW 2 COLUMN 3 
ROW 3 COLUMN 1 
ROW 3 COLUMN 2 
ROW 3 COLUMN 3 
MATRIX-ENTERED 

( 031 ) 


C0S(R[3] ) ; 
SIN ( R[3] ) ; 


0 ; 

-SIN ( [3] ) ; 
COS ( R [3 ] ) ; 
0 ; 

0 ; 

0 ; 

l; 


; cos(r 3 ) 

' -SIN(R 3 ) 
0 


SIN(R 3 ) 

cos(r 3 ) 

0 


(C32) 

ENTERMATRIX(3,3) ; 

ROW 

1 

COLUMN 

1 

C0S(R[2] ) ; 

ROW 

I 

COLUMN 

2 

0; 

ROW 

1 

COLUMN 

3 

-SIN(R[2]) ; 

ROW 

2 

COLUMN 

1 

0; 

ROW 

2 

COLUMN 

2 

l; 

ROW 

2 

COLUMN 

3 

0; 

ROW 

3 

COLUMN 

1 

SIN(R[2] ) ; 

ROW 

3 

COLUMN 

2 

0; 

ROW 

3 

COLUMN 

3 

COS ( R[2 ] ) ; 



MATRIX-ENTERED 


[ COS(Ro) 0 

[ 

(D32 ) [ 0 1 

[ 

[ SIN(R 9 ) 0 

[ 2 

(C33) ENTERMATRIX(3,3) ; 


ROW 

1 

COLUMN 

1 

l; 

ROW 

1 

COLUMN 

2 

0; 

ROW 

1 

COLUMN 

3 

0; 

ROW 

2 

COLUMN 

1 

0; 

ROW 

2 

COLUMN 

2 

COS ( R[ I ] ) ; 

ROW 

2 

COLUMN 

3 

S I N ( R [ 1 ]) ; 

ROW 

3 

COLUMN 

1 

0; 

ROW 

3 

COLUMN 

2 

- S I N ( R [ 1 ] ) ; 

ROW 

3 

COLUMN 

3 

COS ( R [ I ] ) ; 

MATRIX-ENTERED 

(D33) 

[ 

[ 

[ 


[ 

[ 0 

[ 


0 

COS ( R-j ) 
- S I N ( R-j ) 


(C34) ENTERMATRIX(3,1 ) ; 


ROW 1 COLUMN 1 0; 

ROW 2 COLUMN 1 0; 

ROW 3 COLUMN 1 M*G; 

MATRIX-ENTERED 

(D34) 


[ 

[ 

[ 

[ 

[ 


-SIN(R 2 ) 

0 

COS(R 2 ) 


0 ] 

] 

SIN(R 1 ) ] 
1 ] 
COS ( R, ) ] 
1 ] 


0 ] 
] 

0 ] 
] 

, M ] 



The product of these four matrices gives the following column vector of gravitational force com- 
ponents relative to aircraft body axes: 


(C35) (D33) . (D32) . (D31 ). (D34) ; 


(D35) 


[ -SIN(R 9 ) G M ] 
[ 2 ] 
[ SIN(Ri) C0S(R 9 ) G H ] 

[ 1 2 ] 

[ C0S(R-, ) C0S(R 9 ) G M ] 

[ ' 2 ] 


These vector components may be expressed in conventional form by executing the following 
two programming steps, which yield: 

(C36) FOR 1:1 THRU 3 DO R0W[ I ] : FIRST ( ROW ((D35),I))$ 

( C37 ) FOR 1:1 THRU 3 DO ( FG[I] : R0W[I][1 ] , DI SPLAY ( FG[ I ]) )$ 


FG 1 = -SIN(R 2 ) G M 
FG 2 = SIN ( R-j ) C0S(R 2 ) G M 
fg 3 = COS ( R-, ) cos(r 2 ) G M 

where Rj = + 5 Rj), Rf } are equilibrium values, and 8/? ; - are angular perturbations. 


Inertia Forces 

The formulation of the inertia forces involves the determination of the product of an angular 
velocity matrix and a column vector of linear velocity components. This product is the matrix 
equivalent of the familiar vector product cu X V. By adding to the components of this vector, the 
components of linear acceleration relative to aircraft body axes, the components of inertial accelera- 
tion relative to these axes are obtained. The required matrices may be entered and multiplied as 
follows: 


(C38) ENTERMATRIX(3,3) ; 

ROW 1 COLUMN 1 0; 

ROW 1 COLUMN 2 -P[3] ; 

ROW 1 COLUMN 3 P [ 2 ] ; 
ROW 2 COLUMN 1 P[3] ; 
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ROW 2 COLUMN 2 


0; 

ROW 2 COLUMN 3 -P[l]; 

ROW 3 COLUMN 1 - P [ 2 ] ; 

ROW 3 COLUMN 2 P[1 ] ; 

ROW 3 COLUMN 3 0; 


MATRIX-ENTERED 

[ 0 

c 

(D38) [ P 

r j 


(C39) ENTERMATRIX(3,1 ) ; 




] 

] 

] 

] 

] 

] 


ROW 1 COLUMN 1 U[1 ] ; 
ROW 2 COLUMN 1 U[2] ; 
ROW 3 COLUMN 1 U[3] ; 
MATRIX-ENTERED 


(D39 ) 

( C40 ) (D38) . (D39) ; 
( D40 ) 


[ U, ] 
[ 1 1 
[ U ? ] 

[ 2 ] 

[ U^ ] 

[ 3 ] 


£ P 2 U 3 

£ U 1 P 3 
[ P-, u 2 


" U 2 P 3 jj 
" P 1 U 3 j 
- U 1 P 2 } 


( C41 ) FOR 1:1 THRU 3 DO ROW[ I ] : FIRST ( ROW ( ( D40 ) , I ) ) $ 
(C42) FOR 1:1 THRU 3 DO ( C [ I ] : ROW[ I ] [1 DISPLAY (C[ I ] ) )$ 



C 1 P 2 U 3 ‘ U 2 P 3 


C 2 U 1 P 3 " P 1 U 3 


C 3 P ] U 2 U 1 P 2 


A statement of the fact that the z'th component of the linear velocity vector is a function of 
time, requires the use of the DEPENDENCIES function. The use of this function permits the 
system to differentiate the components Uj with respect to time. The remaining two programming 
statements request the system to add the components, multiply the individual sums by the mass M 
of the vehicle, and display the resulting inertial force components FRj as follows: 

(C43) DEPENDENCIES(U(I,T))$ 

(C44) FOR 1:1 THRU 3 DO FR[I]:M*(C[I]+DIFF(U[I],T))$ 

(C45) FOR 1:1 THRU 3 DO D I SPLAY ( FR [ I ] ) $ 

F R, = (1 U, f P 2 U 3 - U 2 P 3 ) M 

FR 2 = U 2 P 1 U 3 + U 1 P 3^ H 

FR 3 * M (" U l P 2 + " U 3 + P 1 U 2 ) 


Resultant Forces 


It only remains to request MACSYMA to combine the aerodynamic, gravitational, and inertia 
forces that were formulated in preceding sections and display the results. The /th component of the 
resultant force will be denoted by FTj where 7/ is the z'th component of thrust. The two program- 
ming steps and the formulated equations follow. 

(C46) FOR 1:1 THRU 3 DO FT[I ] : FR [ I ]- FG [ I ] - FA [ I ]$ 

(C47 ) FOR 1:1 THRU 3 DO DISPLAY(FT[I])$ 


FT, - -(-D F r C0S(A)SIN(B ) + D F r C0S(A)C0S(B) 

i r 2’ L k ] ,L K 


- Dp SIN(A) )DEL(C,,) - (D F ,, C0S(A)SIN(A)SIN 2 (B) 

'3 5 ^K ^ 'o 5U9 


2 2 


- Dr 


,, COS (A)SIN(A)COS(B)SIN(B ) 

> 9 U 1 
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F 1 ,U 2 


F 2’ U 3 


F 3’ U 1 


F 3’ U 3 


F 2’ P 1 


F P 

h ,k 2 


F P 
r 2’ 3 


F P 
h 3 ,r l 


F P 
3’ 3 


F A 
h 2 ,A 1 


F A 
,rt 2 


F A 
r 2’ 3 


F A 
h 3’ M l 


F A 
r 3’ 3 


f 2 ,u 2 


F 3 ,U 1 


F 3 ,U 2 


COS(A)SIN(A)COS(B )SIN(B ) + D p y SIN 2 (A)SIN(B) 

C0S 2 (A)SIN(B) + D F ,, COS(A)SIN(A)COS 2 (B) 
h l’ u l 

SIN 2 (A)C0S(B) + Dp ,, COS 2 (A)COS(B) 
h l’ U 3 

COS(A)SIN(A))DEL(U 3 ) - ( d f 2 ,P 2 COS(A)SIN(A)SIN 2 (B) 
COS(A)SIN(A)COS(B )SIN(B) 

COS(A)SIN(A)COS(B)SIN(B ) + D F p SIN 2 (A)SIN(B) 

r 3’ 2 

■COS 2 (A)SIN(B) + Dp D COS(A)SIN(A)COS 2 (B) 

SIN 2 (A)COS(B) + Dp D C0S 2 (A)C0S(B) 

,P 3 

COS(A)SIN(A})DEL(Po) - (Dp . COS ( A ) S I N ( A ) S I N 2 ( B ) 

J I ^ 5 ^ 

COS(A)SIN(A)COS(B)SIN(B) 


C0S(A)SIN(A)C0S(B )SIN(B ) + D F . SIN 2 (A)SIN(B) 

h 3 ,H 2 

C0S 2 (A)SIN(B) + Dp . C0S(A)SIN(A)C0S 2 (B) 
h l 

SIN 2 (A)C0S(B) + Dp . C0S 2 (A)C0S(B) 

r 3 

C0S(A)SIN(A))DEL(A 3 ) - (" d f 2 ,U-, cos ( a ) sin2 ( B ) 

C0S(A)C0S(B)SIN(B) + D F ,, C0S(A)C0S(B)SIN(B) 

' 1 ,U 1 

SIN(A)SIN(B) + Dp ,, C0S(A)C0S 2 (B) 

' 1 ,U 2 

SIN(A)COS(B))DEL(U 9 ) - (-Dp D C0S(A)SIN 2 (B) 

6 h 2 ,P l 



- Dc D COS(A)COS(B )SIN(B) + D F p COS(A)COS(B )SIN(B ) 

h 2’ P 2 1’ 1 

- Dp D SIN(A)SIN(B) - D F p COS(A)COS 2 (B ) 

h 3’ 1 12 

- Dp D SIN(A)COS(B))DEL(P 9 ) - (-Dp . C0S(A)SIN 2 (B) 

'2^2 c r 2’ M l 

- Dr . COS(A)COS(B ) S I N ( B ) - D F . COS(A)COS(B)SIN(B) 

F 2 ,A 2 r 1 

- D c n SIN(A)SIN(B ) + Dp . COS(A)COS 2 (B) 

P 3’ 1 1’ 2 

- D c n SIN(A)COS(B))DEL(A 9 ) - (Dp .. COS 2 ( A ) S I N 2 ( B ) 

r M 2 ^ ‘ ^ 9 ^ 2 

- Dr ,, C0S 2 (A)C0S(B)SIN(B) - Dp .. C0S 2 (A)C0S(B)SIN(B) 

r l ,u 2 

+ Dp ,, COS (A )SIN(A)SIN(B ) + D p COS(A)SIN(A)SIN(B ) 

5U3 

+ Dr ,, C0S 2 (A)C0S 2 (B) - Dp . C0S(A)SIN(A)C0S(B) 

l "l’ u l 3’ 1 

- Dp [J COS (A)SIN(A)COS(B ) + D p SIN 2 (A) )DEL(U ] ) 

- (Dp D C0S 2 (A)SIN 2 (B) - Dp p C0S 2 (A)C0S(B)SIN(B) 

' 2*^2 ' 2 9 ' 1 

- Dp C0S 2 (A)C0S(B)SIN(B) + d p COS(A)SIN(A)SIN(B) 

P 1’ P 2 3’ 2 

+ Dp D COS (A)SIN(A)SIN(B ) + D F p C0S 2 (A)C0S 2 (B ) 

' 2 ’ ' 3 r 1 

- Dp p C0S(A)SIN(A)C0S(B) - D F p C0S(A)SIN(A)C0S(B) 

h 3 ,p i r 3 


I I 

- Dp p C0S(A)SIN(A)C0S(B ) - D F p C0S(A)SIN(A)C0S(B) 

P 3’ P 1 1’ 3 

+ Dp D SIN 2 (A))DEL(P, ) - (Dp . C0S 2 (A)SIN 2 (B) 

' 3 ,P 3 2’ 2 

- Dp . C0S 2 (A)C0S(B)SIN(B) - Dp . C0S 2 (A)C0S(B)SIN(B) 

+ D F . COS(A)SIN(A)SIN(B ) + D p . COS ( A ) S I N ( A ) S I N ( B ) 

3’ 2 2 ,M 3 
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COS(A)SIN(A)COS(B) 


+ D,- . COS (A)COS (B ) - D P . 

r 3 ’ 1 

- Dp . COS(A)SIN(A)COS(B) + D F . SIN 2 (A) )DEL(A 1 ) 

h l ,m 3 3’ 3 

+ sin(r 2 )gm + (-- u-, + p 2 u 3 - u 2 p 3 )m 

dT 

+ S c COS( A)SIN(B ) - Sp C0S(A)C0S(B ) + S F SIN(A) 

'2 h l 3 

FT 9 = -(Dp r SIN(B) + Dp „ COS(B))DEL(C k ) 
c. ' 1 5 ^ K V 2*^K ^ 

- (-Dp .. SIN(A)SIN 2 (B) - Dp .. S I N ( A ) CO S(B)SIN(B) 

h l ,U 2 2’ 2 

+ D c ,, SIN(A)COS(B)SIN(B) + D F COS(A)SIN(B) 
h l ,U l 1’ 3 


+ D,- .. SIN(A)C0S 2 (B) + Dp .. COS (A)COS(B ) )DEL(lL) 

r ^ ^ U i r 2’ u 3 ^ 

Dp D SIN(A)SIN 2 (B) - Dp p SIN(A)COS(B)SIN(B) 

l" i 5 r 2 ^ 2 9 '2 


- (-D, 


+ Dr 


1 ’ 


+ °F 
- (-D 

+ D 


p SIN(A)COS(B)SIN(B) + D F p C0S(A)SIN(B) 

K 1 h l ,l 3 

p SIN(A)C0S 2 (B) + Dp p cos(a)cos(b))del(pj 
> , h l 2’ 3 

Ip . SIN(A)SIN 2 (B) - Dp . SIN(A)COS(B)SIN(B) 

• i 9 v 2^2 

V n SIN(A)COS(B)SIN(B) + D F „ C0S(A)SIN(B) 

F] ’A] h ] ,A 3 

- . SIN(A)COS 2 (B) + Dp . COS(A)COS(B))DEL(A- 

"2 ,A 1 2’ 3 

Ip .. SIN 2 (B) + Dp .. COS (B)SIN(B) + D F .. C0S(B)SIN(B) 
r-|,U-| r ^ r*| 9^2 

Ip ,, C0S 2 (B))DEL(U 9 ) - (Dp p SIN 2 (B) + Dp p COS(B)SIN(B) 

f 2’ u 2 ' 2 9 ' 1 

+ Dp D C0S(B)SIN(B) + Dp p C0S 2 (B) )DEL(P„) - (D F . SIN 2 (B) 

' 9 ' 2 ' 2 9 ' 2 ^ ' i 


+ D 


(D 


+ D 
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+ D 


+ D 


F 2’ A 1 


F 2 ,A 2 


- D 

- D 

- D 

- D 

- D 

- D 

- D 

- D 

- D 


f 2 ,u 


F r U 


f 2 ,u 

f 2 ,p 


F r p 


f 2 ,p 


f 2 ,a 


F-| ,A 


f 2 ,a 


2 

3 

3 

2 

3 

3 

2 

3 

3 


+ ($- U, 
dT c 


FT 3 = -<- D 


+ D 


F 3’ C K 


- D 


F 2’ U 1 


- D 


F 3’ U 2 


+ D 


F 1 ,U 1 


COS(B)SIN(B) + D P , COS(B )SIN(B ) 

1 ,A 2 

COS 2 (B))DEL(A 2 ) - (-Dp (J cos(a)sin 2 (b) 

COS(A)COS(B)SIN(B) + Dp ,, COS(A)COS(B)SIN(B) 

h r u i 

SIN(A)SIN(B) + Dp .. COS(A)COS 2 (B) 

h 2’ U l 

SIN(A)COS(B ) )DEL(Un ) - (-Dp D COS(A)SIN 2 (B) 

i ' 1 5 ^2 

COS(A)COS(B)SIN(B) + D F D COS(A)COS(B)SIN(B) 

h l 

SIN(A)SIN(B ) + Dp D COS(A)COS 2 (B) 

h 2 ’ 

SIN(A)COS(B))DEL(P 1 ) - (-D p #A COS(A)SIN 2 (B ) 

COS(A)COS(B)SIN(B) + D F . COS(A)COS(B)SIN(B) 

,A 1 

SIN(A)SIN(B) + Dp . COS(A)COS 2 (B) 

h 2 ,A 1 

SIN(A)COS(B) )DEL(A-| ) - SIN(R ] )COS(R 2 )GM 


- P ] U 3 + U-, P 3 )M - S F SIN(B) - Sp^ COS(B) 


F r SIN( A)SIN(B) + Dp r SIN(A)C0S(B ) 
h 2 ,L K l ,b K 

C0S(A))DEL(C k ) - (Dp^ y SIN 2 (A)SIN 2 (B) 

SIN 2 (A)COS(B)SIN(B) - Dp SIN 2 (A)C0S(B)SIN(B) 

,u 2 

COS(A)SIN(A)SIN(B) - D P COS(A)SIN(A)SIN(B) 

' 2 ,u 3 

SIN 2 (A)C0S 2 (B) + Dp 

I '3 ,U 1 


C0S(A)SIN(A)C0S(B) 



+ D p ^ u COS(A)SIN(A)COS(B ) + D p u COS^(A) )DEL(U 3 ) 

13 3 3 

- (D F p SIN 2 (A)SIN 2 (B) - D P D SIN 2 (A)COS(B)SIN(B) 

r 2’ r 2 r 2 » ' i 

- Dp p SI'N 2 (A)COS(B)SIN(B) - Dp D COS(A)SIN(A)SIN(B) 

r l ’ k 2 h 3 ,H 2 

- Dp p COS (A )SIN(A)SIN(B ) + D r D S IN 2 (A )C0S 2 ( B ) 

h 2’ j h l* K l 

+ Dp p COS(A )SIN(A)COS(B ) + D c D C0S(A)SIN(A)C0S 2 (B) 
h 3 ’ K 1 h l’ K l 

- Dp p SIN 2 (A)COS(B) + Dp d C0S 2 (A)C0S(B) 

h l ’ P 3 h 3’ P l 

- Dp p COS(A)SIN(A) )DEL(P, ) - (D P . C0S(A)SIN(A)SIN 2 (B) 

r 3’ j 1 h 2’ rt 2 

- Dp . COS(A)SIN(A)COS(B)SIN(B) 

2 ,M 1 

- Dp . COS (A)SIN(A)COS(B)SIN(B) + D r . SIN 2 (A)SIN(B) 

r l ,a 2 h 2 ,A 3 

- Dp . COS 2 (A)SIN(B) + Dp . COS(A)SIN(A)COS 2 (B) 

3’ 2 

- Dp . SIN 2 (A)C0S(B) + Dp . C0S 2 (A)C0S(B) 

1 ,m 3 h 3’ M l 

- Dp . COS(A)SIN(A))DEL(A n ) + M(-U,P, + -- U, 

r 3 ,M 3 1 1 L dT 6 

+ P ] U 2 ) - COS(R ] )C0S(R 2 )GM + S p SIN(A)SIN(B) 

- Sp SIN(A)C0S(B) - Sp COS(A) 

h 3 

+ Dp D SIN (A ) C0S 2 (B) + Dp p COS(A) COS(B) ) DEL(P ? ) 
h l’ r 2 l "3’ l 2 

- (-Dp . SIN (A) SIN 2 (B) - Dp . SIN (A) COS(B) SIN(B) 

' 2 ‘ 2 9 ™2 

+ Dp . SIN(A) COS(B) SIN(B) + D c . COS(A) SIN(B) 

h 3’ 1 
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+ D F A SIN(A) C0S 2 (B) + D f ^ a ^ COS(A) COS(B)) DEL(A 2 ) 


3’ 2 


- (D,- .. COS(A) SIN(A) SIN 2 (B) - D F ,, COS(A) SIN(A) COS(B) SIN(B) 

*2 > u 2 ' ? 5 u i 


2 1 


- Dr .. COS (A) SIN(A) COS(B) SIN(B) + D F M SIN 2 (A) SIM(B) 

h l’ u 2 2 ,u 3 

- D r „ C0S 2 (A) SIN(B) + D F M COS(A) SIN(A) COS 2 (B) 

h 3’ U 2 h 'l ,U l 

- D F n SIN 2 (A) COS(B) + D F n C0S 2 (A) COS(B) 

l "l ,u 3 3’ 1 

- D r ,, COS(A) SIN (A) ) DEL (U, ) - (D F p COS(A) SIM(A) SIN 2 (B) 

h 3’ 3 2’ 2 

- Dr D COS (A) SIN (A) COS(B) SIN(B) 

' 2 ,h I 

- Dr p COS(A) SIN(A) COS(B) SIN(B) 

h I’ K 2 

+ D F p SIN 2 (A) SIN(B) - D p COS 2 (A) SIN(B) 

h 2’ 3 h 3’ 2 

+ Dr p C0S(A)SIN(A)C0S(B) 

r 3 

+ Dr p C0S 2 (A))DEL(PJ - (Dr A SIN 2 (A)SIN 2 (B) 

' 3 5 ' 3 ^ r 2 5 m 2 


- Dr 


n SIN 2 (A)C0S(B)SIN(B) - Dr A SIN 2 (A)C0S(B)SIN(B) 

I 5 “ 1 r ] ’ H 2 


- Dr A C0S(A)SIN(A)SIN(B) - Dr « COS(A)SIN(A)SIN(B) 

F 3 ,A 2 2’ 3 

+ D r „ SIN 2 (A)C0S 2 (B) + D P C0S(A)SIN(A)C0S(B) 
T fl l 3 ,h 1 

+ D p A C0S(A)SIN(A)C0S(B) + D p A C0S 2 (A) )DEL(A 3 ) 
13 3 3 

- (-Dr ,, SIN(A)SIN 2 (B) - D f u SIN (A)COS(B ) S IN ( B ) 

^2*^1 ' 9 5 u 9 


2 2 


+ D 


r „ SIN(A)C0S(B ) S I N ( B ) + D F ,, C0S(A)SIN(B) 
h l’ u l l ‘3’ 1 


+ d f u SIN(A)C0S 2 (B) + d f u cos(a)cos(b))del(u 2 ) 
12 3 2 



- (-D F p SIN(A)SIN 2 (B) - Dp p SIN(A)COS(B)SIN(B) 
r 2’ 1 2 ,P 2 

+ Dp p SIN(A)COS(B)SIN(B) + D F COS(A)SIN(B) 

1 r 3’ K l 


Special Forms of the Equations of Motion 

In aeronautical studies involving small perturbations about the equilibrium or trim condition, 
the investigator sometimes wants to know how the vehicle will respond if the motion is restricted in 
some way. For example, he might wish to determine vehicle response in the absence of sideslip. 
MACSYMA is well equipped to implement assumptions of this type. By using a substitution com- 
mand, MACSYMA goes through the equations, makes the required substitutions, and displays the 
modified results. For the case of zero sideslip the program requests MACSYMA to make the substi- 
tutions: SIN(B) = 0 and COS(B) = 1 in each force equation. The required substitution and display 
commands and the modified equations assume the following form : 


( C48 ) FOR 1:1 THRU 3 DO FT[I] : SUBST ( [ S I N ( B ) = 0 , COS ( B ) = 1 ] ,FT[I])$ 
(C49 ) FOR 1:1 THRU 3 DO D I SPLAY ( FT [ I ] ) $ 


FT, = -(D c _ COS(A) - Dp r SIN(A))DEL(C„) - (-D F ,, SIN 2 (A) 
I t r L K h 3 ’ L K K h 3 ,U l 


- D 


- .. C0S(A)SIN(A) + Dp C0S(A)SIN(A) 

3 ,U 3 1’ 1 


+ D 


Ip U C0S 2 (A))DEL(U 3 ) - (-d f p SIN 2 (A) - D f 


,p p C0S(A)SIN(A) 
3 ,k 3 

Jp p C0S(A)SIN(A) + Dp p C0S 2 (A))DEL(P-) - (-D F . SIN 2 (A) 

r 3 3’ i 

- Dp „ C0S(A)SIN(A) 

'^ 5^3 


+ D c 


+ 


1 


. C0S(A)SIN(A) 

’ i 


+ D 


- (D 

- D 


l F A C0S 2 (A))DEL(A 3 ) - (Dp^ 5lJ COS(A) - D p ^ ^ SIN ( A ) )DEL (U 2 ) 

ip p COS (A) - D f 
1 ,k 2 


v P SIN(A) )DEL(P ? ) - (Dp . COS(A) 

r 2> r 2 ^ r i » “2 

Ip n SIN(A))DEL(A 9 ) - (Dp ,, SIN 2 (A) - Dp C0S(A)SIN(A) 
' 3 ,m 2 c 3’ u 3 h 3’ u l 


Dp C0S(A)SIN(A) 
h l ,U 3 


+ D 


■p U C0S 2 ( A) ) DEL ( U-j ) - (D f 
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- D 


+ D 


c p COS(A)SIN(A) - D F p COS(A)SIN(A) 

1-3 , h 1 r r K 3 

P p COS 2 (A))DEL(P 1 ) - (Dp . SIN 2 (A) - Dp . COS(A)SIN(A) 
h l* H l 1 h 3’ A 3 h 3 ,H l 

- D. . COS(A)SIN(A) + Dp . C0S 2 (A) )DEL(A, ) + SIN(R ? )GM 

h l ,A 3 h l’ A l 

+ (-- U 1 + P 2 U 3 - U 2 P 3 )M + S p SIN(A) - S F COS(A) 
dT 3 1 

FT 2 = - d f 2 ,c k d *< c k> - (D F 2 ,u 1 sin < a > + D F 2 ,u 3 C0S(A))DEL(U 3 ) 

- (Dp p SIN(A) + Dp p C0S(A))DEL(P3) - (D p . SIN(A) 

F 2’ P 1 F 2’ P 3 2 ,M 1 

+ D F 2 ,A 3 COS(A))DEL(A 3 ) - DEL(U 2 ) - d F 2 > p 2 del(p 2 ) 

- Dp . DEL(Ap) - (Dp ,, COS(A) - D p ,, SIN(A) )DEL(U-. ) 

h 2 ’ A 2 ^ r 2 ,u l r 2 ,u 3 


- (D 


F 2’ P 1 


COS(A) 


- D p p SIN(A) ) DEL ( P-j ) - (D F ^ A ^ 


COS(A) 


- Dp A SIN(A))DEL(A 1 ) - SIN(R ] )COS(R 2 )GM 
2 3 

+ ( 5 - u 2 - p,u 3 + u,p 3 )m - s 

dT c 


FT- = -(Dp r SIN(A) + Dp r COS(A) ) DEL ( C ^ ) 

o r -| 5 1 ^ 


3 K 


- (Dp .. SIN 2 (A) + Dp .. COS(A)SIN(A) + D F COS(A)SIN(A) 

F 1’ U 1 h 3’ U l 1’ 3 

+ Dp ,, COS 2 (A))DEL(U,) - (Dp p SIN 2 (A) + D p p COS(A)SIN(A) 
F 3’ u 3 3 F 1’ P 1 h 3’ P l 

+ Dp D COS(A)SIN(A) + Dp D C0S 2 (A) )DEL(P J - (D F fl SIN 2 (A) 

F 1 ’ P 3 ' 3 ,P 3 3 1’ 1 

+ Dp . COS(A)SIN (A) + Dp . C0S(A)SIN(A) 

h 3’ A l 1’ 3 

+ Dp A C0S 2 (A))DEL(A 3 ) - (Dp ^ SIN(A) + D p ^ COS(A) )DEL(U 2 ) 
3 3 1 2 3 2 



- (D F jP SIN(A) + D F ^ C0S(A))DEL(P 2 ) - (D p ^ SIN(A) 

1 2 3 2 1 2 

+ Dp A COS (A) )DEL( A ? ) - (-D r SIN 2 (A) - D P COS(A)SIN(A) 
r 3’ H 2 * h l’ u 3 h 3 ,U 3 

+ Dp^^ COS(A)SIN(A) + D p ^ COS 2 (A) )DEL(U ] ) 

- (-Dp SIN 2 (A) - D F p COS(A)SIN(A) + D P D COS(A)SIN(A) 

r l ,h 3 h 3’ 3 h 1’ F l 

+ Dp COS 2 (A))DEL(P-. ) - (-Dp . SIN 2 (A) 
h 3 ,K l 1 h l’ A 3 


D r 


A COS(A)SIN(A) + Dp . COS (A)SIN(A) 
!’ m 3 h r M i 


+ Dp . C0S 2 (A))DEL(A 1 ) + M(-U,P, + -- U, + P,U„ 
h 3’ A l 1 1 6 dT J 1 L 


- C0S(R )COS(R )GM - S P SIN(A) - S F COS(A) 

' * r l h 3 

In addition to the zero sideslip condition, the investigator might wish to determine vehicle 
response when the angle of attack A is limited to small values. For this condition the program 
would request MACSYMA to make the substitution S1N(A) = A. Moreover, if the angle of attack 
were sufficiently small, the program would request MACSYMA to make the additional substitution 
COS(A) = 1 . 

In this case, the required substitutions and display commands give rise to the following modi- 
fied equations: 


(C50) FOR 1:1 THRU 3 DO FT[I ] : SUBST( [SIN(A )=A,C0S( A)=l ] ,FT[ I] )$ 
(C51 ) FOR 1:1 THRU 3 DO DISPLAY(FT[I])$ 


FT 1 (Di V c |< " ° F 3’ C k A)DEL(C K ) ' (_D F 3 ,U 1 A " D F 3 ,U 3 


A 


+ D 
+ D 


F 1 ,U 1 A + °F 1 ,U 3 )DEL(U 3 ) " (_D F 3 ,P 1 " "F Q ,P,, " ' “Ft.P 


3 1 

A 2 - D r 


A + D F n A 
3”3 h l ,h l 


F 1 , P 3 )DEL(F>3) ‘ ( ' D F 3 ,A 1 A " D F 3 ,A 3 A + D F 1 ,A 1 A 


l’°l 


+ D F 1 ,A 3 ) DEL ^ 3 ) - (°F r U 2 - D F 3 ,U 2 A)DEL(U 2 ) - (D^ 
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V p A)DEL(P ? ) - (D f - D f A)DEL(A ? ) - (D p ,, 
r ^ s r ^ *- r i 5 m ^ v 3 9 


Dp y A - D 

i- 3 , u ] 


F r u 3 A + D FrU 1 )DEL(u i ) ' (d f 3 ,p 


A - D f p A 
3 '3*^1 


- Dp p A + Dp p ) DEL ( P ) - (Dp „ A^ - Dp . A - Dp . A 

1’ 3 r 1 3’ 3 3’ 1 1’ 3 

+ Dp . )DEL(A-, ) + SIN(R,)GM + (-- U, + P 9 U, - UJ>)M + S c A - Sp 
1 c dT ^ 6 c 6 h 3 1 

FT 2 = - D F 2 ,C K DEL(C k) - ( D F 25Ui A + D F 2 ,U 3 ) DEL(U 3 ) - (D F 2 , Pl A 

+ D F 2 ,P 3 > DEL < P 3> - (°F 2 ,A 1 A + - D F 2 ,U 2 DEL(U 2> 

“ D F 2 ,P 2 DEL ^ P 2^ " D F 2 ,A 2 DEL ^ A 2^ " ^ D F 2 ,U 1 ” D F 2 ,U 3 a ) del ( u i) 

" (DE 2’ P l" DF 2’ P 3 A)DEL(Pl) " (Dp 2’ A l " ° F 2’ A 3 A)DEL(Al) 

- SIN(R ] )C0S(R 2 )GM + (^- U 2 - P 1 U 3 + U ] P 3 )M - S p ^ 

FT 3--« D F 1 .C |C A + D F 3 .C K > DEL(C K>-< D F 1 .U 1 a2 + D F 3 ,U 1 A 

+ D F 1 ,U 3 A + D F 3 ,U 3 ) [)E1 -( U 3> - 'V, ft2 + Vi fl + V 3 A 
+ D F 3 > P 3 ) DEL < P 3 ) - (“FpA, + “Fj.A, A + D F r A 3 A 
+ D F 3 ,A 3 > DEL(A 3) - (D F r U 2 A + D F 3 .U 2 > DEL ' U 2) - (D F,,P 2 fl 

+ D F 3 ,P 2 > DEL < P 2) - < D F r A 2 A + D F 3 ,A 2 ) DEL ( A 2> - <- D F r U 3 A2 

-Dp ,, A + Dp A + Dp .. )DEL(U 1 ) - (-D F p A 2 - D p p A 
F 3 ,U 3 F r U l F 3 ,U 1 1 F l’ p 3 h 3V 

+ D. A+D p ) DEL ( P-i ) - (-Dp A A 2 - D A + D A 

F 1’ P 1 F 3’ P 1 1 F 1’ A 3 ' 3’ A 3 V A 1 
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+ D F A )DEL(A 1 ) + M(-U-, P 2 + " U 3 + P ] U 2 ) “ COS(R 1 )C0S(R 2 )GM 

- S F A - S P 
h l h 3 

Examination of these equations reveals the existence of terms such as A 2 . If it is assumed that 
second-order terms in A are negligible, a program statement instructing MACSYMA to make the 
substitution A 2 = 0 would simplify the equations as follows: 

(C52) FOR 1:1 THRU 3 DO FT[I ] : SUBST( [A**2=0] ,FT[I] )$ 

(C53) FOR 1:1 THRU 3 DO DISPLAY(FT[I])$ 


FT l = -< d f,,c k - D F 3 ,C K a > del < c k> - <- d f 3 ,u 3 a + °F 1 ,u 1 a 


+ D F 3 ,U 3 > DEL ^ U 3 ) - <- D F 3 ,P 3 A + °F r P, A + D F 1 ,P 3 ’ DEL ( P 3> 


- (-D P . A + Dp . A + Dp . )DEL(A,) - (D F 

r 3 5 3 r -j 9 r\ -j * "1 9 *0 J i 5 


u. 


- D 


A - D 


F 1 ,A 3 


D F 3 ,U 2 A ) DEL(U 2) - (D F r P 2 - D F 3 ,P 2 A)DEL(P 2) - ( V 
F 3 ,A 2 A )° EL ( A 2 ) - (- D F 3 ,U 1 A - D F r U 3 A + D F rUl ) DEU 

“ ^ _d f 3 ,p 1 a “ D F v P 3 a + °F 1 ,P n ) DEL ( P 1 ) - (- d f 3 ,a 1 

+ Dr- . ) DEL ( A-, ) + SIN(R 9 )GM + (-- U-, + P 9 U 9 - U 9 P 9 )M+S F A - S P 
h i ’ A i 1 ^ dT J 3 1 

FT 2 = - D F 2 ,C k DEL ^ C k) - ( D F 2 ,U 1 A + D F 2 ,U 3 ) DEL(U 3) - (D F 2 ,P-, A 

+ D F 2 ,P 3 ) DEL ( P 3) " ( D F 2 , Ai A + D F 2 ,A 3 > DEL < A 3> - D F 2 ,U 2 DEL(U 2) 

- d f 2 ,p 2 DEL(p 2 ) - d f 2 ,a 2 del (V - - d f 2 ,u 3 a > del < u i> 

- ( d f 2 , Pi - d f 2 ,p 3 A > DEL < p i> - (°f 2 ,a 1 “ d f 2 ,a 3 a)del(a i } 

- SIN(R 1 )C0S(R 2 )GM + U 2 - P 1 U 3 + U 1 P 3 )M - S p ^ 
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FT 3 " "^ D F-,,C k A + D F 3 , C|< ) DEL ( C K> " ^ D F 3 ,U- ] A + D F p U 3 A + D F 3 ,U, 

DEL(Uo) - (D F p A + D f p A + D f p )DEL(P-) - (D p - A 
j r 3’ r l r l’ r 3 r 3’ 3 3’ 1 

+ Vs A + D f 3 ,a 3 ^ del ^ A 3^ - (D F r u 2 A + d f 3 ,u 2 ) del ( u 2 ) 

- (Dp p A + Dp p )DEL(P ? ) - (Dp . A + Dp . )DEL(A ? ) 

• 1 >» 2 V 2 1’ 2 r 3 9/ *2 ^ 


^" Df 3’ U 3 A + DF 1 ,U 1 A + ^' Df 3’ P 


A + D f p A 
3 h ,f l 


+ Dp p ) DEL ( P-, ) - (-Dp A A + Dp fl A + D p . )DEL(A ] ) 

h 3’ 1 3’ 3 1’ 1 3’ 1 

+ M(-U-|P 2 + -- U 3 + P-jUg) - C0S(R-| )C0S(R 2 )GM - S F A - S F 
dT r l r 3 

Additional simplifications are possible if it is assumed that angular velocity perturbations are 
negligible. This assumption can be implemented by again using the substitution command, which 
yields the following greatly simplified equations: 

(C54) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
FT [ I ] : SUBST( [DEL ( P [ J ] ) =0] , FT [ I ] ) $ 

(C56) FOR 1:1 THRU 3 DO DISPLAY(FT[I])$ 

FT 1 = -< d f v c k - d f 3 ,c k a > del(c k> - <- d f 3 ,u 3 a + “f,^ fl 

+ Dp )DEL(U,) - { -Dp A+Dp . A + D )DEL(A ) 

r 1 ’ u 3 ^ r 3 ,m 3 r ] "I ,n 3 ^ 


- (° F u - D 

h l’ U 2 


F ^ A)DEL(U 2 ) - (D F jA - Dp A)DEL(A 2 ) 
3 2 1 2 3 2 


- ("D c I, A - Dp .. A + D 
3’ 1 1’ 3 


F 0 ,u n " "f 15 u q '■ ’ "f 1 ,u 1 )del(u i ) " (_d f 3 ,a 1 " ^F 15 A 


A - Dr- A 
1 ,rt 3 


+ Dp A ) DEL ( A-, ) + SIN(R 2 )GM + (jj T U ] + P 2 U 3 - U 2 P 3 )M+S p A -S p 
i i 3 i 


FT 2 _D F 2 ,C k DEL ( C k^ ^ D F 2 ,U 1 A + D F 2 ,U 3 ^ DEL ^ U 3' j " ^°F 2 ,A 1 A 

+ D f „ )DEL(AJ - Dp (| DEL(U 2 ) - D p fl DEL(A 2 ) 

^ 2 5 "'S ^ 2’ 2 ^ * 2 5 ^*2 
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(D 


f 2>Ui " d f 2 ,u 3 a )° el ( u i) ' ^ d f 2 ,a 1 - d f 2 ,a 3 A )° EL ( A i 


SIN(R-, )COS(R 2 )GM + (-- U 2 - P-|U 3 + U-,P 3 )M - S p 

dT 2 


n 3 “ + D F 3 .C K > DEL(C K> - '“Fj.U, 


A + D f u A + D F u ) 
r l ,u 3 r 3’ u 3 


q i ~ v^c a n ■ u p a A + Dp n )DEL(A 3 ) - (Dp 
o r ^,m-| r l ,M 3 1 3’ n 3 ° ' 1 9 


DEL(U^) - (D c . A + D 
3’ i 


U, 


+ D F 3 ,U 2 ^ DEL ^ U 2^ ^ Dfr v A 2 A + D F 3 ,A 2 ^ DEL ^ A 2^ ( “ D F 3 ,U 3 A 

+ [ Vi a + d f 3 ,u 1 ) del < u iJ “ ( -°f 3 ,a ; 

+ Dp A JDEKA^ + M(-U 1 P 2 + ~ u 3 + P 1 U 2 ) 


A + Dp n A 
'3 1’ 1 


- COS(R-j )cos(r 2 )gm - Sp 


A - S r 


Finally, it may be of interest to consider the effect of omitting the linear acceleration terms. 
By comparing the response of the system with and without acceleration perturbations, the influence 
of these perturbations can be determined. Again, a simple substitution command is all that is 
required to implement the assumption that DEL(A|)=0. Execution of this command yields the 
modified equations as follows; 


, FOR 1 : 1 IH'RU 3 DO FOR j:l THRU 3 DO 

FT f I ] : SUB ST ([DEL(A[J]j = Q] ) FT[I])S 


CSS) FOR in THRU 3 DO DISPLAY ( FTC I ])S 


H 1 = °F,,C„ A ) del C c k) "l ,,i.i M ■ “F, ,U, 


1 - K 


A + ", ,, A 

3 ,w 3 h l- 1 


# D F „ )DEL(U,: 
r l -3 


< d f, ,u. 


:, F„U, MOEMUj) - (-d f>iUi a 

C, 


3 -1 


= Dp y A + Dp_ _ y_ ) DEL (U 1 ) * SIN(R 2 )GM + (2- ^ + PgUg - "J 1 ,)M 


1 V A “ s F 

3 r l 
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FT 2 = - D F 2 ,C k DEL(C K> ' < D F 2 ,U, A + D F 2 ,U 3 > DEL < U 3> - D F 2 ,U 2 DEL(U 2> 

- (D f LJ - D F u A ) DEL ( U-j ) - SIN(R 1 )COS(R 2 )GM 
2 1 2 3 

+ <aT U 2 - P 1 U 3 + U 1 P 3> M - S F, 


ft 3 “ " <d f 1 ,c k fl + d f 3 ,c k ,del(c k ) ' (d f 3 ,u 1 A + D F,,U, A + D F,,U. 1 


1 3 


3 3 


DEL < U 3> - ( D F r U 2 fl + D F 3 ,U 2 > DEL < U 25 - (- D F 3 ,U 3 A + D F,,l 

+ D F 3 ,U 1 )DEL(U 1 ) + H( - U 1 P 2 + “ U 3 + P 1 U 2> 

- COS(R, )COS(R 9 )GM - S F A - S P 
1 L h 3 


A 


Thrust Forces 


It should be noted that the thrust forces FT[ appearing on the left-hand side of these equations 
are the resultant of a number of thrust generating systems, each contributing a thrust vector T n . 
Each thrust vector is referred to a thrust axes system X n l with origin at the point of application of 
the thrust vector. The axes are chosen such that each thrust vector coincides with the X n x axis of 
the system. Moreover, each thrust vector is then transformed to a coordinate system Y n l which has 
the same origin as the thrust axes, but is parallel to the body axes system. Finally, the components 
of thrust in the Y n ‘ system of axes are transformed to the body axes system, which has its origin at 
the center of gravity of the aircraft. Each thrust axis X n l is related to the Y n i system by the follow- 
ing transformation equations (see sketch): 



Y n l =X n ' cos(K n ) cos(P n ) 

Y n = X n cos( ^« } Sin(/ V ’ 

V=-V s w«) 


Hence, the components of the thrust vector T n in the Yj system of coordinates are 


(13) 


9 V 


T 

n ’ 


dY r 


9 V 


T ■ 
n ’ 


9 V 


(14) 
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These are also the components of thrust in the y i system of coordinates, which has its origin at 
the center of gravity of the aircraft. The thrust components due to all thrust generating systems are 
obtained by summing the right-hand side of the following equation 


<_ a v 


T = 


9 V 


(15) 


The expanded form of equation (15), when summed over n will yield the resultant thrust com- 
ponents. When the number of thrust generating systems is known, the components V can be formu- 
lated and displayed by using equation (15), and executing the following two commands, which 
yield the components contributed by the nth thrust generating system. These are 


(Cl ) Y[1 ,N]:X[1 ,N]*C0S(K[N])*C0S(P[N])$ 

(C2) Y[2 ,N] : X[1 ,N]*COS(K[N])*SIN(P[N])$ 

(C3) Y[3,N]:-X[1,N]*SIN(K[N])$ 

(C4) FOR I THRU 3 DO T[I]:DIFF(Y[I,N],X[1 ,N],1 )*T[N]$ 
(C5) FOR 1:1 THRU 3 DO DISPLAY (T[I])$ 


T, = T n COS ( K n ) C0S(P N ) \ 

T 2 = T n C0S(K n ) SIN(P n ) > (16) 

T 3 ' -T„ SIN(K n ) J 

Determination of the Geographical Location of Aircraft 

In order to determine the geographical location of an aircraft relative to some initial location, 
it is necessary to transform the components of the aircraft’s velocity vector from aircraft body axes 
to a system of earth-fixed axes. The transformed components can then be integrated to find the 
location of the aircraft as a function of time. The product of the three rotation matrices (D33), 
(D32), and (D31), which were used to transform the gravity vector from an earth-fixed axes system 
to aircraft body axes, may be transposed and used to transform the aircraft velocity components 
to an earth-fixed system. If the column vector (D39) of aircraft velocity components is premulti- 
plied by the transposed matrix, the velocity components relative to the earth-fixed system are 
obtained as follows: 

(C56) TRANSPOSE( ( D33) . (D32) . (D31 ) ). (D39) ; 
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( D56) MATRIX([U 3 (SIN(R 1 )SIN(R 3 ) + COS(R-j )SIN(R 2 )COS(R 3 ) ) 

+ U 2 (SIN(R 1 )SIN(R 2 )C0S(R 3 ) - COS(R 1 )SIN(R 3 )) + u 1 cos(r 2 )cos(r 3 )], 

[U 2 (SIN(R 1 )SIN(R 2 )SIN(R 3 ) + C0S(R 1 )C0S(R 3 )) 

+ U 3 (COS(R 1 )SIN(R 2 )SIN(R 3 ) - SIN(R 1 )COS(R 3 )) + U n COS(R 2 )SIN(R 3 )], 

[-U 1 SIN(R 2 ) + U 2 SIN(R 1 )COS(R 2 ) + U 3 C0S(R 1 )C0S(R 2 )]) 

If the components X[ relative to the Earth-fixed system be denoted by DX(, execution of the 
following programming steps will ensure that the required velocity components are displayed in 
conventional form. 

(C57 ) FOR 1:1 THRU 3 DO ROW [ I ] : FI RST ( ROW ( ( D56 ) , I ) )$ 

(C58) FOR 1:1 THRU 3 DO ( DX [ I ] : ROW [ I ] [ 1 ]. DISPLAY ( DX[ I ] ) ) ; 

DX-j = U 3 (SIN(R-|)SIN(R 3 ) + C0S(R-! )SIN(R 2 )C0S(R 3 )) 

+ U 2 (SIN(R 1 )SIN(P 2 )COS(R 3 ) - C0S(R 1 )SIN(R 3 )) + U 1 C0S(R 2 )C0S(R 3 ) 

DX 2 = U 2 (SIN(R 1 )SIN(R 2 )SIN(R 3 ) + C0S(R 1 )C0S(R 3 ) ) 

+ U 3 (C0S (R-, )SIN(R 2 )SIN(R 3 ) - SIN(R 1 )C0S(R 3 )) + U 1 COS(R 2 )SIN(R 3 ) 

DX 3 = -U-J SIN( R 2 ) + U 2 SIN(R 1 )C0S(R 2 ) + U 3 C0S(R-j )cos(r 2 ) 

Integration of these velocity components will yield the required coordinates of the aircraft 
relative to a set of earth-fixed reference axes. These are 

X E Uxi EO + J DXmdt 

where X^q are the initial values of the coordinates in the earth-fixed reference frame. 

MOMENTS 


Transformation of Static Moments 

The static aerodynamic moments obey the same transformation law as the static aerodynamic 
forces; that is, if Sm h denotes a static moment in the X frame of reference, and SMj denotes the 
corresponding transformed moment in the Y reference frame, then 


(17) 


SM t = 


dY' 

dX n 



where Y = Y(X) is obtained from the displayed output (D6) and reentered here to facilitate the 
formulation of the moment equations. Given the transformation equations (D6), the transformed 
aerodynamic static moments are obtained by expanding equation (17). The three programming 
steps used to transform the static forces may again be employed to transform the static moments. 
The simple program and the displayed results are 


(Cl ) Y[1]:X[1]*C0S{A)*C0S(B)-X[2]*C0S(A)*SIN(B)-X[3]*SIN(A)$ 

(C2) Y[2]:X[1]*SIN(B)+X[2]*C0S(B)$ 

(C3) Y[3]:X[1]*SIN(A)*C0S(B)-X[2]*SIN(A)*SIN(B)+X[3]*C0S(A)$ 

(C4) SM[ I ] : =0$ 

(C5) FOR 1:1 THRU 3 DO FOR N:1 THRU 3 DO 
SM[I]:SM[I]+DIFF(Y[I],X[N])*S[M[N]]$ 

(C6) FOR 1:1 THRU 3 DO DISPLAY (SM[I])$ 

SM-j = -S^ C0S(A)SIN(B) + S M C0S(A)C0S(B) - S M SIN(A) 

SM 2 = S M SIN(B) + S M COS(B) 

SM 3 - -S^ SIN(A)SIN(B) + SIN(A)C0S(B) + S M COS(A) 

Transformation of Control Moment Derivatives 

The control moment derivatives obey the same transformation law as the static moments; that 
is, if Dm h ,Ck denotes the nth control moment derivative with respect to the ATh control surface as 
measured in the X reference frame, and TDj q denotes the corresponding transformed derivative in 
the Y frame, then 


TD 


I,C 


dY l 

dx n 


D 


' M rv C k 


(18) 


where Y = Y(X) is again obtained from the displayed output (D6). 

As in the preceding section, the transformed control derivatives are obtained by expanding the 
transformation law (18) given the transformation equations (D6). The transformed derivatives are 
obtained by executing the following simple program, which has exactly the same form as the pro- 
gram used to transform the static moments. These are 
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(C7) TD[I,C]:=0$ 

(C8) FOR 1:1 THRU 3 DO FOR N:1 THRU 3 DO 
TD[I,C]:TD[I,C]+DIFF(Y[I],X[N])*D[M[N],C[K]]$ 

(C9) FOR 1:1 THRU 3 DO DISPLAY(TD[I ,C])$ 


TD 


1,C 


TD 


3,C 



TD 



COS(A )SIN (B ) + D m r COS(A)COS(B ) - D M r 

K V L K n 3 ,L 

2.C = y C|< SIN ^ + D M 2 ,C k C0S ^ 

SIN(A)SIN(B) + D m SIN(A)COS(B ) + D^ 


SIN(A) 


COS(A) 


The corresponding control moments are obtained by multiplying the control derivatives by the 
appropriate control increments DEL{Cj^). The following two programming steps are sufficient to 
formulate the required moments. These are denoted by CM j in the displayed output. 


(CIO) FOR 1:1 THRU 3 DO CM[I] :TD[ I ,C]*DEL(C[K] )$ 
(Cl 1 ) FOR 1:1 THRU 3 DO DISPLAY(CM[I ])$ 


CM-, 


(-D m r C0S(A)SIN(B) + D M r C0S(A)C0S(B) 
2 ,l K m 1 ,l K 


- D m q SIN(A))DEL(C k ) 
3 ’ K 


CM. 


CM. 


(D m jC sin (B ) + D m c C0S(B))DEL(C k ) 

1 ’ K 2 ’ K 

(-D m r SIN(A)SIN(B) + D m r SIN(A)C0S(B) 
2 ,l K V L K 


+ D m ?c C0S(A))DEL(C k ) 
3 K 


Moments Produced by Linear Velocity Perturbations 

The next step in the formulation involves the determination of the aerodynamic moments pro- 
duced when an aircraft is subjected to linear velocity perturbations DEL (Uj). Before these moments 
can be determined, the aerodynamic stability derivatives with respect to linear velocity components 
must be transformed from wind or wind-tunnel stability axes to body axes. For a detailed discus- 
sion of the transformation of these derivatives, the reader is referred to equations (4) through (10). 
The program used for the transformation of force derivatives can be used in this case also. In this 
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application, the aerodynamic stability derivative of the / th moment with respect to the /th velocity 
component will be denoted by Z?M/, Ui - The corresponding transformed derivatives are denoted by 
TDMi,Uj- When the program is rewritten to accommodate the notational changes required for this 
application, it assumes the following form: 

(Cl 2) TDU [ I ,J]:=0$ 

(Cl 3) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
FOR R: 1 THRU 3 DO FOR N:1 THRU 3 DO 

TDU [I , J ] :TDU[I , J]+DIFF( Y[I ] ,X[R])*DIFF(Y[J],X[N])*D[M[R],U[N]]$ 

It only remains to multiply the transformed derivatives by the appropriate velocity increments 
to obtain the required moments, which are denoted by MDU /. The next three programming steps 
instruct MACSYMA to evaluate and display the moments produced by linear velocity perturbations. 
These are 


(Cl 4) MDU [ I ] : =0$ 

(Cl 5) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
MDU[I]:MDU[I]+TDU[I,J]*DEL(U[J])$ 

(Cl 6) FOR 1:1 THRU 3 DO D I S PL A Y ( MDU [ I ] ) $ 


MDU -j = (D m u C0S(A)SIN(A)SIN (B) 

- D m u C0S(A)SIN(A)C0S(B)SIN(B) 

- D m .. C0S(A)SIN(A)C0S(B )SIN(B) + D M ,, SIN 2 (A)SIN(b; 

ri-1 ,1)2 3*^2 

- D m u C0S 2 (A)SIN(B) + D M C0S(A)SIN(A)C0S 2 (B) 

- D m u SIN 2 (A)C0S(B) + d Mi ,u 3 C0S 2 (A)C0S(B) 


- D m y C0S(A)SIN(A))DEL(U 3 ) + (-D m ^u^ COS(A)SIN^B) 


- D k 


u C0S(A)C0S(B )SIN(B ) + D^ y C0S(A)C0S(B)SIN(B) 


- D m ^ 5U ^ SIN(A)SIN(B) + D M C0S(A)C0S^(B) 

- D m ^ ;U ^ SIN(A)C0S(B))DEL(U 2 ) + (D M25 u 2 cos 2 (a)sin 2 (b) 
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- COS 2 (A)COS(B)SIN(B) - D M JJ C0S 2 (A)C0S(B)SIN(B) 

+ D m ^ COS(A)SIN(A)SIN(B) + D m jL) COS(A)SIN(A)SIN(B) 

3 2 2 5 3 

+ COS 2 (A)COS 2 (B) - COS(A)SIN(A)COS(B) 

- D m ^ COS(A)SIN(A)COS(B) + D m jLJ SIN 2 (A))DEL(U 1 ) 

13 3’ 3 

MDU 2 = (-D^^ SIN(A)SIN 2 (B) - D M y S I N ( A ) COS (B)SIN(B) 

+ SIN(A)COS(B)SIN(B) + D Mi#u cos(a)sin(b) 

+ D m SIN(A)COS 2 (B) + D COS(A)COS(B))DEL(U,) 

2 5 "| 2 ? ^3 ^ 

+ (D M 1 ,U 1 SIn2 ( B ) + COS(B)SIN(B) + D„ COS(B)SIN(B) 

+ C0S 2 (B))DEL(U 2 ) + (-D M y C0S(A)SIN 2 (B) 

- D m ^ COS (A )COS (B )SIN(B ) + ^ COS(A)COS(B)SIN(B) 

2 2 1 ] 

- D m ?u SIN(A)SIN(B) + D m ^ C0S(A)C0S 2 (B) 

1 3 2 5 ] 

- % \j SIN(A)COS(B))DEL(U 1 ) 

MDU 3 = (D M 2 ,U 2 si N 2 (A)SIN 2 (B) - D M SIN 2 (A)COS(B)SIN(B) 

- D m u SIN 2 (A)COS(B)SIN(B) - D m COS(A)SIN(A)SIM(B) 

1 ’ 2 3 ,u 2 

- D M ^ COS(A)SIN(A)SIN(B ) + D M ;U SIN 2 (A)C0S 2 (B) 

2 3 1 ’ 1 

+ D m u COS(A)SIN(A)COS(B) + D m COS(A)SIM(A)COS(B) 

3’ 1 1 ’ u 3 

+ D m ^ COS 2 (A))DEL(U 3 ) + (-D m y SIN(A)SIN 2 (B) 

3 3 2 ’ 1 



I 


- D^ ^ SIN(A)COS(B)SIN(B) + D m y SIN(A)COS(B)SIN(B) 

+ D M ^ cos (a )sin(b ) + D M ^ SIN(A)COS 2 (B) 

3 1 12 

+ D m ^ C0S(A)C0S(B))DEL(U 2 ) + (D m y cos(a)sin(a)sin 2 (b) 
3 2 2 2 

- D M y COS(A)SIN(A)COS(B)SIN(B) 

- D m ^u 2 COS(A)SIN(A)COS(B)SIN(B) + D m u sin 2 (a)sin(b) 

- D m u cos 2 (a)sin(b) + D m u C0S(A)SIN(A)C0S 2 (B) 

3 2 11 

- D M ^ SIN 2 (A)C0S(B) + 0 C0S 2 (A)C0S(B) 

- D M y COS(A)SIN(A))DEL(U 1 ) 


Moments Produced by Angular Velocity Perturbations 

The program used in the preceding section can, with suitable notational changes, be used to 
formulate the moments produced by angular velocity perturbations. However, whereas in the pre- 
ceding application the required moments were obtained by multiplying the transformed aerody- 
namic stability derivatives by linear velocity increments, in the present case the transformed deriva- 
tives must be multiplied by angular velocity increments. In view of these similarities, the following 
program and displayed moments will be presented without further comment, except to point out 
that the aerodynamic stability derivatives of the / th moment with respect to the / th angular velocity 
component are denoted by DfylpPj- The corresponding transformed derivatives are denoted by 
TDMi'Pj, and the resulting moments by MDPj. 


(Cl 7 ) TDP [ I ,J]:=0$ 

(C18) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
FOR R: 1 THRU 3 DO FOR N:1 THRU 3 DO 

TDP[ I , J] : TDP[I , J]+DIFF( Y[I] ,X[R] )*DIFF(Y[J ],X[N])*D[M[R],P[N]$ 
(Cl 9) MDP [ I ] : =0$ 

(C21 ) FOR 1:1 THRU 3 DO FOR J : 1 THRU 3 DO 
MDP[ I] : MDP[I ]+TDP[I , J ]*DEL(P[J ] )$ 

(C22) FOR 1:1 THRU 3 DO DISPLAY(MDP[I])$ 
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MDP-, = (D m p COS(A)SIN(A)SIN 2 (B) 

- D m D COS(A)SIN(A)COS(B)SIN(B) 

n 2’ K l 

- D m p COS(A)SIN(A)COS(B)SIN(B) + D M p SIN 2 (A )SIN(B ) 

12 3 2 

- D M p^ COS 2 (A)SIN(B) + p C0S(A)SIN(A)C0S 2 (B) 

- D M p SIN 2 (A)C0S(B) + D M p C0S 2 (A)C0S(B) 

- D M p C0S(A)SIN(A))DEL(P 3 ) + (-D M p COS(A)SIN 2 (B) 

- D m COS(A)COS(B)SIN(B) + D m d COS(A)COS(B)SIN(B) 

^ 2 9 ' 2 1 5 ' 1 

- D M p SIN(A)SIN(B) + D M p COS (A )C0S 2 (B ) 

- D M P 2 SIN(A)COS(B))DEL(P 2 ) + (D M p cos 2 (a)sin 2 (b) 

- p^ COS 2 (A)COS(B)SIN(B) - p^ COS 2 (A)COS(B)SIN(B) 


+ D M p^ COS(A)SIN(A)SIN(B) + D m p COS(A)SIN(A)SIN(B) 

+ d m p cos 2 (a)cos 2 (b) - D m d COS(A)SIN(A)COS(B) 

"v 1 3’ 1 


- D m 5 p COS(A)SIN(A)COS(B) + d m 5 p SIIn (A) )DEL(P ] ) 
13 3 3 


MDP 2 = (-D m p SIN(A)SIN (B) - D m SIN(A)COS(B)SIN(B) 


+ D m SIN(A)COS(B)SIN(B) + D p COS(A)SIN(B ) 

"r i V 3 

+ D M p SIN(A)C0S 2 (B) + D M p C0S(A)C0S(B))DEL(P 3 ) 
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+ (D M p SIN 2 (B) + D M p COS(B)SIN(B) + D M p COS(B)SIN(B) 
+ d m p C0S 2 (B))DEL(P 2 ) + (-d m p^ C0S(A)SIN 2 (B) 

- D m p COS(A)COS(B)SIN(B) + D m COS(A)COS(B)SIN(B) 

- D M p SIN(A)SIN(B) + D^p COS(A)COS 2 (B) 

- D M p SIN(A)COS(B))DEL(P 1 ) 

MDP 3 = (D m p SIN 2 (A)SIN 2 (B) - D m p sin 2 (a)cos(b)sin(b) 

- D m p SIN 2 (A)C0S(B)SIN(B) - p^ COS(A)SIN(A)SIN(B) 

- D m p COS(A)SIN(A)SIN(B) + D m p sin 2 (a)cos 2 (b) 

+ D m p COS(A)SIN(A)COS(B) + p COS(A)SIN(A)COS(B) 

+ D M p COS 2 (A))DEL(P 3 ) + (-d m p SIN(A)SIN 2 (B) 

- D m p SIN(A)COS(B)SIN(B) + D m p^ SIN(A)COS(B)SIN(B) 

+ D M 5 p COS(A)SIN(B ) + D M >p sin(a)cos 2 (b) 

3 1 12 

+ D m 5 p COS(A)COS(B))DEL(P 2 ) + (D m jP COS(A)SIN(A)SIN 2 (B) 

3 2 2 5 2 

- d m p COS(A)SIN(A)COS(B)SIN(B) 

2 ’ K 1 

- D m p COS(A)SIN(A)COS(B)SIN(B ) + d m p SIN 2 (A)SIN(B) 

1 ’ 2 2 ’ r 3 



- D M D COS 2 (A)SIN(B) + D m p C0S(A)SIN(A)C0S 2 (B) 

m 3 ,k 2 1 ’ 1 

- D M p SIN 2 (A)COS(B) + D M p COS 2 (A)COS(B) 

- D M p COS(A)SIN(A))DEL(P 1 ) 


The same procedure may be used to formulate the aerodynamic moments produced by linear 
and angular accelerations. These moments will not be included here, since the cases considered so 
far are sufficient to demonstrate the facility with which symbolic mathematical computation can be 
used to formulate and transform aerodynamic moments. 

Inertia Moments 

The formulation of inertia moments involves the determination of the product of an angular 
velocity matrix, a matrix of inertia coefficients, and a column vector of angular velocity compo- 
nents. This product is the matrix equivalent of the familiar vector product cu X //, where is the 
angular velocity vector and H is the angular momentum vector. By adding to the components of this 
vector, a vector which represents the rate of change of angular momentum relative to the moving 
body axes, the inertial moments relative to these axes are obtained. The rate of change of angular 
momentum relative to the moving body axes may be expressed as the product of the inertia matrix 
and a column vector of angular acceleration components. The required matrices may be entered and 
multiplied as follows: The first matrix to be entered is the inertia matrix, with elements J { j. It is 
entered by typing the statement ENTERMATRIX(3.3) and responding to the system’s request for 
elements. 

(C23) ENTERMATRIX( 3,3) ; 

ROW 1 COLUMN 1 J[1 ,1 ]; 

ROW 1 COLUMN 2 J[1 ,2]; 

ROW 1 COLUMN 3 J[1 , 3] ; 

ROW 2 COLUMN 1 J[2,l]; 

ROW 2 COLUMN 2 J[2,2]; 

ROW 2 COLUMN 3 J[2,3]; 

ROW 3 COLUMN 1 J[3,l]; 

ROW 3 COLUMN 2 J[3,2]; 


48 



ROW 3 COLUMN 3 J [ 3 , 3 ] ; 


MATRIX-ENTERED 

: J i. 

I 

J l. 

2 

T 

(D23) 

■ °2, 

I 

J 2, 

2 

J 2, 


■ 

1 

J 3, 

2 

J 3, 


A statement of the fact that the z'th component of the angular velocity vector is a function of 
time requires the use of the DEPENDENCIES function. The use of this function permits the system 
to differentiate the components Pj with respect to time, and to enter the resulting acceleration 
components in the form of a column vector as follows: 

(C24 ) DEPENDENCIES(P(I,T))$ 

(C25) ENTERMATRIX(3,1 ) ; 

ROW 1 COLUMN 1 DIFF(P[1 ],T) ; 

ROW 2 COLUMN 1 DIFF(P[2] ,T) ; 

ROW 3 COLUMN 1 D IFF ( P [3] »T ) ; 


MATRIX-ENTERED 


(025) 


dT 


dT 


dT 


P 


P 


P 


] 

1 ] 
] 
] 
] 

2 ] 
] 
] 
] 

3 ] 

] 


The angular velocity matrix and a column vector of angular velocity components are entered 

next 


(C26 ) ENTERMATRIX(3,3) ; 

ROW 1 COLUMN 1 0; 

ROW 1 COLUMN 2 -P[3] ; 
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ROW 1 COLUMN 3 


P[2]; 

ROW 2 COLUMN 1 P[3]; 

ROW 2 COLUMN 2 0; 

ROW 2 COLUMN 3 -P[l]; 

ROW 3 COLUMN 1 -P[2]; 

ROW 3 COLUMN 2 P[1 ] ; 

ROW 3 COLUMN 3 0; 

MATRIX-ENTERED 

( D26 ) 


(C27 ) ENTERMATRIX(3, 1 ) ; 



‘ ] 
- p ’ 3 

o ] 
] 


ROW 

1 

COLUMN 

1 

PD]; 

ROW 

2 

COLUMN 

1 

P [2] ; 

ROW 

3 

COLUMN 

1 

P[3] ; 

MATRIX 

-ENTERED 



(D27) 


P H 

D 

D 


These four matrices are now combined to yield a column vector of inertia moments relative to 
aircraft body axes. 

(C28) ( (D23) . (D25)+(D26) . (D23) . (D27 ) ) ; 


(D28) MATRIX([J (^- P ) + J (^- P ) + J (4_ P ) 
i,J d T J '>* dT ^ 1,1 dT 1 



[J 2,3 P 3 } + J 2,2 P 2 ) + J 2,l P 1 ) ~ P 1 (P 3 J 3,3 + 


P 0 J 


dT 


dT 


2 3,2 


+ P 1 J 3,1 ) + P 1 J 3,1 ) + P 3 (J 1,3 P 3 + J 1,2 P 2 + P 1 J 1,1^’ ^3,3 *V 
+ J 3,2 P 2 ^ + J 3,l P l } + P 1^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 


dT 


dT 


P 2^ J 1 ,3 P 3 + J 1 ,2 P 2 + P 1 J 1 ,1 ^ 

The next two programming steps enable the system to express these inertia moments in con- 
ventional functional form. 

(C29) FOR 1:1 THRU 3 DO ROW [I ] : FIRST (ROW ( (D28 ) , I ) )$ 

(C30) FOR 1:1 THRU 3 DO ( IM[I ] : ROW [I ][1 ] , DISPLAY ( IM[I ] ) )$ 

IM 1 = J 1 , 3 P 3^ + J l,2 P 2 ^ + J l,l P l^ + P 2^ P 3 J 3,3 

+ 3 5 2 + P 1 J 3,1^ " P 3^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 

IM 2 = J 2,3 <" P 3> + J 2,2 (" P 2> + J 2,l <" P 1 > " P 1< P 3 J 3,3 
+ ? 2 J 3 t Z + P ! J 3 5 ]) + P 3^ J 1,3 P 3 + J 1,2 P 2 + P 1 J 1,1^ 

™3 = J 3,3 P 3^ + J 3,2 (~ d ~ P 2^ + J 3,l V + P 1^ J 2,3 P 3 


+ P 2^2 , 2 + P 1 J 2,1^ P 2^ J 1,3 P 3 + J 1,2 P 2 + P 1 J 1,1^ 


Resultant Moments 

It only remains to request MACSYMA to combine the aerodynamic and inertia moments 
which have been formulated in preceding sections and to display the results. The z'th component of 
the resultant moment will be denoted by TM t , where TM[ is the ith component of the moment due 
to thrust. The two programming steps and the formulated equations follow. 

(C31 ) FOR 1:1 THRU 3 DO TM[I]: IM[I]-SM[I]-CM[I]-MDU[I]-MDP[I]$ 

(C32) FOR 1:1 THRU 3 DO DISPLAY(TM[I])$ 
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TM 


1 = -(-D m c COS(A)SIN(B) + D M jC COS(A)COS(B) 

2 5 K I s K 

D m jC SIN(A))DEL(C k ) - (D m sU COS(A)SIN(A)SIN 2 (B) 

3 ’ K 2 ’ 2 

D m u COS(A)SIN(A)COS(B)SIN(B) 

D m .. COS(A)SIN(A)COS(B)SIN(B) + D M .. SIN 2 (A)SIN(B) 

Im-j 5 u 2 ^ 2 9 ^2 

D m u COS 2 (A)SIN(B) + D m u C0S(A)SIN(A)C0S 2 (B) 

u SIN 2 (A)COS(B ) + D M y C0S 2 (A)C0S(B) 

D m u eOS(A)SIN(A))DEL(U 3 ) - (D M p C0S(A)SIN(A)SIN 2 (B) 

D m d COS(A)SIN(A)COS(B)SIN(B) 

2’ P 1 

COS (A)SIN(A') CGS (B)SIN(B) + SIN 2 (A)SIN(B ) 

D m p C0S 2 (A)SIN(B) + D m p C0S(A)SIN(A)C0S 2 (B) 

2 3 11 

D M p SIN 2 (A)C0S(B) + D M p C0S 2 (A)C0S(B) 

d m 3 ,p 3 COS(A)SIN(A))DEL(P 3 ) - (-D M u C0S(A)SIN 2 (B) 

D m C0S(A)C0S(B)SIN(B) + D m .. C0S(A)C0S(B)SIN(B ) 
n 2 5 u 2 *^ 1 5 ^1 

D M y SIN(A)SIN(B) + D M y C0S(A)C0S 2 (B) 

D m ^ SIN(A)C0S(B))DEL(U 2 ) - (-D m ?p C0S(A)SIN 2 (B) 

3 2 2 1 

C0S(A)C0S(B)SIN(B) + D m p C0S(A)C0S(B)SIN(B) 
D m ^p^ SIN(A)SIN(B) + D M P 2 C0S(A)C0S 2 (B) 



- d m p SIN(A)COS(B))DEL(P, 
n 3’ 2 £ 


) - (d m u C0S 2 (A)SIN 2 (B) 

) M u C0S 2 (A)C0S(B)SIN(B) 

l M ,, COS(A)SIN(A)SIN(B ) 
n 2’ u 3 


- D m .. COS 2 (A)COS(B)SIN(B) - D ™'' 2 

** 2 9 ^ i 

+ D M u COS(A)SIN(A)SIN(B) + D 
+ D M u C0S 2 (A)C0S 2 (B) - D m y oOS(A)SIN(A)COS(B) 

- D m COS(A)SIN(A)COS(B) + d m ,, SIN 2 (A))DEL(U,) 

T 3 M 3* 3 1 

- (D M p 2 COS 2 (A)SIN 2 (B) - D M p COS 2 (A )COS(B )SIN(B ) 

- D M p COS 2 (A)COS(B)SIN(B) + D m p COS(A)SIN(A)SIN(B) 

+ p COS(A)SIN(A)SIN(B) + D m p COS 2 (A)COS 2 (B) 

- D m p COS(A)SIN(A)COS(B) - D m d COS(A)SIN(A)COS(B) 

3’^i M r K 3 

+ D M p SIN 2 (A))DEL(P 1 ) + S M COS (A)SIN(B ) - S M COS(A)COS(B) 

+ S SIN(A) + J (^- P ) + J (^- P ) + J (^- P ) 

"3 1 ’ J dT J 1 dT ^ 1 ’ 1 dT 1 

+ *V P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1^ " P 3^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2 , 1 ^ 

™2 = _(D M ,C SIN(B) + °M ,c cos(b))del(c k ) 

1 ’ K 2 ’ K 

- (-D^ u ^ SIN(A)SIN 2 (B) - D m lj SIN(A)COS(B)SIN(B) 

+ °M 1 ,U 1 SIN(A)COS(B)SIN(B) + D m u COS(A)SIN(B) 

+ sin(a)cos 2 (b) + d m y C0S(A)C0S(B))DEL(U 3 ) 

- (- D M 1? P 2 SIN(A)SIN 2 (B ) - D M p SIN(A)COS(B)SIN(B) 
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+ D m d SIN(A)COS(B)SIN(B) + d m p COS(A)SIN(B) 

1 ’ 3 


+ D M jP sin(a)cos 2 (b) + D m >p COS(A)COS(B))DEL(P 3 ) 

2 1 2 3 

- (D M sin 2 (b) + COS(B)SIN(B) + COS(B)SIN(B) 

+ D M y C0S 2 (B))DEL(U 2 ) - sin 2 (b) + COS(B)SIN(B) 

+ D m p COS(B)SIN(B) + D M ^p^ cos 2 (b))del(p 2 ) 

- (-D M U2 cos (a)sin 2 ( B ) - COS(A)COS(B)SIN(B) 

+ D m u COS(A)COS(B)SIN(B) - D m ^ SIN(A)SIN(B) 

11 13 

+ D m ^ COS(A)COS 2 (B) - D m sU SIN(A)COS(B))DEL(U 1 ) 

2 1 2 3 

- (-D m p 2 COS(A)SIN 2 (B) - D m ^p^ COS(A)COS(B)SIN(B) 


+ D m p COS(A)COS(B)SIN(B) 
1 ,f l 


) - D M p SIN(A)SIN(B) 

"l ,k 3 

+ D M 5 p COS(A)COS 2 (B) - D m 5 p SIN(A)COS(B ) )DEL(P 1 ) - s m sir 

2 1 2 3 1 

- S„ 2 COS(B) + J 2i3 (^- P 3 > + J 2,2<" P 2> +J 2,l<~ P 1 > 

" P 1^ P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1^ + l V J l,3 P 3 + J 1,2 P 2 + P 1 J 1,1^ 

™3 = _(_D M C SIN ( A ) SIN ( B ) + d m ,c SIN ( A ) C0S ( B ) 

2 ’ K 1 ’ K 

+ D m c COS(A))DEL(C k ) - (d m sU SIN 2 (A)SIN 2 (B) 

3 ’ K 2 ’ 2 

- D m .. SIN 2 (A)COS(B)SIN(B) - d m sin 2 (a)cos(b)sin(b) 

1^2 5 u -j * 5 ^2 

- D m u COS(A)SIN(A)SIN(B) - D m u cos(a)sin(a)sin(b) 


) - 


2 1 ( -- P l } 

L ’ 1 HT 1 
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COS(A)SIN(A)COS(B ) 


+ d m,,u, SIN 2 (A)C0S 2 (B) + a„^ 

+ D m .. COS(A)SIN(A)COS(B) + D M C0S 2 (A) )DEL(U 3 ) 

M p u 3 3 ,u 3 

- (D M D SIN 2 (A)SIN 2 (B) - D m p SIN 2 (A)COS(B)SIN(B) 

MgsTg 2’ r l 

- D m d SIN 2 (A)COS(B)SIN(B) - d m p COS(A)SIN(A)SIN(B) 

m 1» P 2 m 3* 2 

- D m p COS(A)SIN(A)SIN(B) + d m p SIN 2 (A)COS 2 (B) 

M 2’ P 3 1’ 1 

+ D M p COS(A)SIN(A)COS(B) + D M p COS(A)SIN(A)COS(B) 

+ D M p C0S 2 (A))DEL(P 3 ) - (-D^^ SIN(A)SIN 2 (B) 

- D m (J SIN(A)COS(B)SIN(B) + D M y SIN(A)COS(B)SIN(B) 

+ d m y COS(A)SIN(B ) + D M SIN(A)COS 2 (B) 

+ d m y COS(A)COS(B))DEL(U 2 ) - (-d m p SIN(A)SIN 2 (B) 

- D m p SIN(A)COS(B)SIN(B) + p SIN(A)COS(B)SIN(B) 

+ D M p COS(A)SIN(B) + D M p^ SIN(A)C0S 2 (B) 

+ D M p^ C0S(A)C0S(B))DEL(P 2 ) - (D m y COS(A)SIN(A)SIN 2 (B) 

- D M COS(A)SIN(A)COS(B ) S I N ( B ) 

- y^ COS(A)SIN(A)COS(B)SIN(B) + y SIN 2 (A)SIN(B) 

- d m y COS 2 (A)SIN(B) + d m y C0S(A)SIN(A)C0S 2 (B) 

- D M y SIN 2 (A)C0S(B ) + D M y COS 2 (A)COS(B ) 



- D m u COS(A)SIN(A))DEL(U 1 ) - C0S(A)SIN(A)SIN 2 (B) 


- D m d COS(A)SIN(A)COS(B)SIN(B) 

“Vl 


- D m d COS(A)SIN(A)COS(B)SIN(B) + D m p SIN 2 (A)SIN(B) 

m 1’ p 2 n 2* 3 

- d m d cos 2 (a)sin(b) + d m p C0S(A)SIN(A)C0S 2 (B) 

M 3 ,p 2 1’ 1 

- D M p sin 2 (a)cos(b) + D M 5 p C0S 2 (A)C0S(B) 

13 3 1 

- D M COS(A)SIN(A))DEL(P 1 ) + SIN(A)SIN(B) 

- S M SIN(A)COS (B ) - S M COS(A) + J- J-- PJ + J. J-- P ? ) 

M 3 ’ dT J dT *- 

+ J 3,l^‘ V + P 1^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 

" P 2^ J 1,3 P 3 + J 1,2 P 2 + P 1 J 1,1 ) 


Special Forms of the Moment Equation 

As in the case of the force equations, the investigator sometimes wishes to modify the moment 
equations to determine how the vehicle will respond if the motion is restricted in some way. For the 
case of zero sideslip, MACSYMA goes through the equations, makes the appropriate substitutions, 
and displays the modified results. The zero sideslip condition requires that SIN(B) = 0 and 
COS(B) = 1. The substitution and display statements required to implement this assumption and 
the modified moment equations assume the following form: 

(C33) FOR 1:1 THRU 3 DO 
TM[I] :SUBST ( [SIN(B)=0, COS ( B ) =1 ] ,TM[I] )$ 

(C34) FOR 1:1 THRU 3 DO DISPLAY(TM[I])$ 


,c K cos(fl) " D M 3 ,C k SIN(A))DEL(C k ) 

- (-D h y SIN Z (A) - D m _y COS (A)SIN(A) + D H _y COS ( A )S I N ( A) 

3 1 3 3 11 

+ D M ;U C0S 2 (A))DEL(U 3 ) - (-D m ;P sin 2 (a) - D M 5 p C0S(A)SIN{A) 
1 5 3 3 5 1 33 
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+ D m }P COS(A)SIN(A) + D m jP C0S^(A))DEL(P 3 ) - (D m sU cos(a) 
11 13 12 

- D M 3 ,u 2 £: '-(a»del(u 2 ) - (D M)1 P 2 COS(fl) - SIN(A))DEU(P 2 ) 

- (D m sU sin 2 (a) - D m >(j COS(A)SIN(A) - D h >l| COS(A)SIN(A) 

3 3 3 1 1 3 

+ D m ^ C0S 2 (A))DEL(U 1 ) - (D m jP SIN 2 (A) - D m COS(A)SIN(A) 

1 1 3 3 3 1 

- D M jP COS(A)SIN(A) + D m jP C0S 2 (A))DEL(P 1 ) + S M SIN(A) 

13 11 3 

- S COS(A) + J (L P ) + J p ) + J (d- p ) 

1 l,J dT J 1 dT 6 1 » 1 dT 1 

+ P 2 ( p 3 j 3,3 + P 2°3,2 + P 1 J 3,1^ " P 3^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 

™2 = - d m 2 ,c k del < c k> - < d m 2 ,u, SIN(fl) + D H 2 ,U 3 COS(A))DEL(U 3 ) 

- < D M 2 ,P, SIN < A > + d M 2 ,P 3 COS(A))DEL(P 3 ) - DEL(U 2 ) 

- D M 2 ,P 2 DEL < P 2> - ( £) M 2 ,U 1 C0S < A > - D M 2 ,U 3 SIN(A))DEL(U,) 

- < d m 2 ,p, cos < a > - d h 2 ,p 3 SIN(A))DEL(P,) ♦ J 2>3 (1 P 3 ) 

+ J 2,2<~ p 2 ) - S H 2 + P ,> - p ,( p 3 j 3,3 + P 2 J 3,2 + P 1 J 3,1> 

+ P 3 (J 1 ,3 P 3 + °1 ,2 P 2 + P 1 J 1 ,1 * 

™3 = -<Vc„ SIN<A) + D M,.C„ COS(A))DEL(C k ) - (D SIN 2 (A) 

I IN J K 11 

+ D m sU COS(A)SIN(A) + D m sU COS(A)SIN(A) 

3 1 1 5 3 

+ D m }U C0S 2 (A))DEL(U 3 ) - (D m 5 p sin 2 (a) + D m jP COS(A)SIN(A) 
33 1 5 1 3 9 1 

+ D M jP cos (A)SIN(A) + D m 5 p C0S 2 (A))DEL(P 3 ) - (D m jLJ SIN(A) 
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+ D m u COS(A))DEL(U 2 ) - (D m 5 p SIN(A) + D M 5 p cos(a))del(p 2 ) 

3 2 1 2 3 2 

■ (_D M r U 3 SIn2 ( A ) - D M 3 ,U 3 COS(A)SIN(A) + COS(A)SIN(A) 

+ D M ^ C0S 2 (A))DEL(U 1 ) - (-d m 5 p sin 2 (a) - D m 5 p COS(A)SIN(A) 

3 1 13 3 3 

+ D M p COS (A)SIN(A) + D M p COS 2 (A))DEL(P 1 ) - SIN(A) 

- S M COS (A) + j 3 p 3 ) + p 2 ) + j 3 P,) 

M 3 dT J dT ’ dT 

+ P 1^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2 , 1 ^ " P 2^ J 1 ,3 P 3 + J 1,2 P 2 + P 1 J 1,1^ 

In addition to the zero sideslip condition, the investigator frequently wishes to determine 
vehicle response when the angle of attack is limited to small values. For this condition MACSYMA 
would implement the assumption that SIN(A) = A. Moreover, if the angle of attack were suffi- 
ciently small, the program would request MACSYMA to make the additional substitution 
COS(A) = 1. 

In this case, the required substitution and display statements give rise to the following modi- 
fied moment equations: 

(C35) FOR 1:1 THRU 3 DO TM[ I ] : SUBST( [SIN ( A )=A ,C0S ( A ) =1 ] ,TM[I])$ 

(C36 ) FOR 1:1 THRU 3 DO D I S PLA Y ( TM [ I ] ) $ 


r c K “ d m 3 ,c k “ ^" d m 3 ,u 1 a " d m 3 ,u 3 a 


+ D M 1 ,U 1 A + V 3 > DEL ( U 3) ^ ° M 3’ P 1 A ° m d A + D “ n A 


'M 3 ,P 3 " ' “M-, ,P 1 


'1 51 1 


+ D Mi 5 P 3 )D E L( p 3 ) - (D^- D M35U2 A)DEL(U 2 ) - (D^ 
- D M 3 ,P 2 a)de L (p 2 ) - (D„ 3>U3 a 2 - d„ 3;Ui a - A 


+ D M 1 ,U 1 )DEL(U 1 ) " (D M 3 ,P, n ■ U M,, P n n - u M lt P 


A P A ” P A 

3 m 3 ’ 1 "T 3 


+ D )DEL(P ) + S A + J (4- P ) + J (4- P ) + J (4- P ) 
M r r l 3 dT J 1 dT ^ 1 ’ 1 dT 1 
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S M + P 2^ P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1^ *V J 2,3 P 3 + P 2 J 2,2 


+ P 1 J 2,1 ) 

™2 = " D M 2 » C k DEL(C K } " ^ D M 2 ,U 1 A + D M 2 ,U 3 ^ DEL ^ U 3^ " ^ D M 2 ,P-, A 
+ D M 2 ,P 3 > DEL ( P 3) - D M 2 ,U 2 DEL ( U 2> " D M 2 ,P 2 DEL < P 2> - ( D M 2jUi 
- D H 2 ,U 3 A > DEL (’ J 1 ) - - D M 2 ,P 3 A > DEL <P 1 ) + J 2,3<“ P 3> 

+ J 2,2^j P 2* " S H 2 + J 2,d"" P U ' P A P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,P 


+ P 3^ J 1 ,3 P 3 + J 1 ,2 P 2 + P 1 J 1 ,1 ^ 


M 3’ U 1 


™3 = - ( VC K A+ D M 3 ,C k ) DEL ^ C k) - ( D M 1?U] A + D 
+ D M r U 3 + D M3 ,U 3 ) DEL < U 3) - < D H 1 ,P 1 ft2 + Vl A + D H r P ; 


+ D M 3 ,P 3 ) DEL < P 3> - <Vu 2 A + D H 3 ,P 2 > DEL ( U 2> - <V 2 A 


+ D H 3 ,P 2 ) DEL ( P 2> - <-V, 


1 ’3 


A U A + U 

3 ,u 3 M l’ 1 


A 


+ D m ,, ) DEL ( U i ) - ( -D 


A - D„ n A + D, 


'M 3 ,U 1 / “-"'l^ v ^.P," “M^.P, ■’ ' “M n ,P 


3 ’ 3 


n 


+ D M 3 ,P 1 )DEE(P l ) - S H, A + J 3 , 3 <^ P 3 > - S M, + J 3 , 2 ( :; P 2 > 


dT 


+ J 3j( HT P l^ + V J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 


dT 

P 2 ^ J 1 , 3 P 3 + J 1 , 2 P 2 + P 1 J 1 ,1 ^ 

Examination of these equations reveals the existence of terms such as A 2 . If it is assumed that 
second-order terms in A are negligible, a program statement instructing MACSYMA to make the 
substitution A 2 = 0 would simplify the moment equations as follows: 

(C 37 ) FOR 1:1 THRU 3 DO 
TM[ I ] : SUBST ( [A** 2 = 0 ] ,TM[I ] )$ 


(C 39 ) FOR 1:1 THRU 3 DO DI SPLAY (TM[ I ] ) $ 
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™1 = • ( 0 M 1 ,C 


K - D H 3 ,C k A > DEL < C K> - (- D M 3 ,U 3 ft + “m^U, fl 


+ D Ml ,U 3 > DEL < U 3> - ( ' D M 3 ,P 3 fl + D M,, P 1 A + Vp 3 > DEL(P 3> 

- ( D M r U 2 - D H 3 ,U 2 A > DEL < U 2 ) - ( D M r P 2 - D M 3 ,P 2 A > DEL < P 2> 

■ ( ' D M 3 ,U 1 A " D M 1 ,U 3 A + %.U 1 > DEL < U 1 > " ( ' D M 3 ,P 1 A ' D M r P 3 A 
+ D H 1 ,P 1 ) DEL ( P 1» + S M 3 A+J 1,3(" P 3> +J 1,2<^ P 2> 

+ P l ! " S M 1 + P 2 (P 3 J 3,3 + P 2°3,2 + P l J 3,l' 

- p 3 ( J 2j3 p 3 + P 2 J 2,2 + P 1 J 2,l ' 


™2 = -°M 2 ,C k DEL < C K> - < D M 2 ,U, A + V,U,) DEL(U 3> - < %,?, 


+ D H 2 ,P 3 > DEL ( P 3> - D M 2 ,U 2 DEL(U 2> - D M 2> P 2 DEL(P 2> - <V, 

- D H 2 ,U 3 A > DEL <V - <V, - D H 2 ,P 3 A > DEE < P 1> + J 2,3 ( ~ P 3> 

+ J 2.2 ( 3 j P 2* ' S M 2 + °2.1 ( 3 y P l ) " P 1 (P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1 ) 
+ P 3^ J 1 ,3 P 3 + J 1 ,2 P 2 + P 1 J 1 ,1 ^ 


™3 "^ D M 15 C k A + D M 3 ,C k ^ DEL ^ C K^ A + D M-,,U 3 A 

+ - % P] A + D M r P 3 A + D M 3 ,P 3 ) DEL ( P 3 ) 

^ D M 1? U 2 A + D M 3 ,U 2 ^ DEL ^ U 2^ " ^ D M r P 2 A + D M 3 ,P 2 ^ DEL ^ P 2^ 


1 d m 3 ,u 3 a + d m 1 ,u 1 a + d m 3 ,u 1 )del ^ u i ) ■ ^ _d m 3 ,p 3 a + d m 1 ,p 1 a 


+ D M 3 ,P 1 ) DEL ( P l> - S H, ft + J 3,3<~ P 3> - S M 3 + J 3,2<~ P 2> 
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■> + *V + P 1 ^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1^ 

” P 2^ J 1 ,3 P 3 + J 1 ,2 P 2 + P 1 J 1 ,1 ^ 

An additional simplification is possible if the assumption that angular velocity perturbations 
are negligible is a valid one. Implementation of the assumption that DEL{Pj) = 0 yields the follow- 
ing greatly simplified equations: 

(C40) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
TM[I ] : SUBST { [DEL(P[J])=0] ,TM [ I ] ) $ 

(C41 ) FOR 1:1 THRU 3 DO DISPLAY(TM[I])$ 


™1 ( D M r C 


K " Dm 3’ C K A)DEL(C K ) " ^ ° M 3’ U 3 A + V U 1 A 


+ D M n ,U 3 )DEL(U 3> ' \,U 2 - \,U 2 A > DEL <V - <- D M 3 , Ul A 


V u 3 A 


rd 


+ + S M 3 A + J l,3 ( d ~ P 3 y ' u l,2 v :; r 2 


) + J n P 9 ) 
dT 


+ Jl ’^dT ^ S + P 2^ P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1^ 
P 3^ J 2,3 P 3 + P 2 J 2,2 + P 1 J 2,1 ^ 


TM„ = -D M . DEL ( C w ) - (D 


'2 " m 2 ,C K vu M 9 ,U-, A + D M 0 ,U 0 ^ DEL ^ U 3^ 


2 ’1 


2 3 


- D M 2 ,U 2 DEL(U 2) - - D M 2 ,U 3 A > DEL (V + J 2,3 ( at P 3) 

+ J 2,2^" P 2^ " S M 2 + J 2,l^~ V " P 1^ P 3 J 3,3 + P 2 J 3,2 + P 1 J 3,1^ 
+ P 3^ J 1 ,3 P 3 + J 1 ,2 P 2 + P 1 J 1 ,1 ^ 


™3 "^ D M 1 ,C k A + D M 3 ,C k ) DEL ( C k) " ^°M 3 ,U 1 
+ D M 3 ,U 3 ^ DEL ^ U 3^ " ( D M-,,U 2 A + D M 3 ,U 2 ^ DEL ^ U 2 


A + D, 


M 1 ,U 3 


) - (-D 


M 3’ U 3 


+ D, 


q.u, A + D M 3 ,U,> DEL ( U 1> - S M, A + J 3,3<“ P 3 


) - s. 
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Thrust Moments 

As indicated previously, the thrust moments TMj appearing on the left-hand side of these equa- 
tions are the resultant of the moments produced by a number of thrust generating systems. Equa- 
tions (13) relate the thrust axes coordinates X n 1 to the coordinate system Y n ‘, which has the same 
origin as the thrust axes but is parallel to the body axes systems. 

To facilitate the formulation, equations (13) are entered here. 

(Cl ) Y[1 ,N]:X[1 ,N]*C0S( K[N] )*C0S (P[N] )$ 

(C2 ) Y[2 ,N] : X[1 ,N]*C0S( K[N] )*SIN( P[N] )$ 

(C3 ) Y[3,N]:-X[1 , N]*SIN( K[N] )$ 

The point of application of the /Ah thrust vector relative to the body axes system, with origin 
at the center of gravity, has components (L l n ,L 2 ,«.T 3 n ). The components of the nth thrust vector 
in this coordinate system are given by equations (14). The product of the position matrix with ele- 
ments (Zq ' U ,L 2 n ,L 3 n ) and a column vector of thrust components can be processed as follows. 

First enter the (3.3) position matrix, element by element, as requested by MACSYMA. Next 
enter the (3.1) column vector of thrust components in the same manner. When the matrices are 
entered, the displayed form of each matrix assumes the conventional textbook form 

(C4) ENTERMATRIX(3,3) ; 

ROW 1 COLUMN 1 0; 

ROW 1 COLUMN 2 -L[3,N]; 

ROW 1 COLUMN 3 L[2,N]; 

ROW 2 COLUMN 1 L[3,N] ; 

ROW 2 COLUMN 2 0; 

ROW 2 COLUMN 3 -L[l , N] ; 

ROW 3 COLUMN 1 -L[2,N]; 

ROW 3 COLUMN 2 L[1,N]; 

ROW 3 COLUMN 3 0; 
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MATRIX- ENTERED 



[ o 
C 

-L 3,N 

(04) 

[ L 3,N 

0 


[ “ L 2,N 

h,N 


( C 5 ) ENTERMATRIX(3,1 ) ; 


ROW 1 COLUMN 1 T[N]*DIFF(Y[1,N],X[1 , N]) ; 

ROW 2 COLUMN 1 T[N]*DIFF( Y[2,N],X[1 ,N] ) ; 

ROW 3 COLUMN 1 T[N]*DIFF(Y[3,N] ,X[1 , N] ) ; 

MATRIX-ENTERED 

[ t n cos(k n ) cos(p n ) ] 

(05) [ T n COS(K n ) SIN(P n ) ] 

[ -t n sin(k n ) ] 

By requesting the system to multiply these two matrices, the following product matrix is 
obtained: 


(C6) (D4).(D5); 


(06) 


; -l 3>n t n cos(k n ) sin(p n ) - l 2jN t n SIN(K n ) ] 

: l 3iN T n C0S(K n ) C0S(P n ) ♦ L 1jN T n SIN(K n ) ] 

; h,N T N C0S < K N> SIN(P N> " L 2,N T N C0S<K N> C0 S(P n ) ] 


In order to express this column vector of thrust moments in conventional functional form, the 
following two programming steps are required: 

(C7 ) FOR 1:1 THRU 3 DO R0W[I] : FIRST( R0W( (D6) , I ) ) ; 

(C8) FOR 1:1 THRU 3 DO (TM[I]:R0W[I][1],DISPLAY(TM[I]))$ 

TM-| = -L^^ COS(K^) SIN(P^j) - ^2^ T^ SIN(K^) 

™2 = L 3, N T N C0S ( K N } C0S(F V + L 1,N T N SIN(K N } 

™3 = L 1,N T N C0S( V SIN(P N } " L 2,N T N cos <S } C0S(P N ) 
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These equations give the moments produced by the nth thrust vector. When the number of 
thrust generating systems is known, these equations can be summed on n to obtain the total thrust 
moments. 


Spatial Orientation in Terms of the Direction Cosines 

The differential equations for the direction cosines can be obtained by first entering a (3,1) 
column vector of direction cosines, with elements Dj x , D / 2 , and D/ 3 , where / can assume the values 
1 ,2,3, and by premultiplying this vector by the angular velocity matrix. This operation ^equivalent 
to the vector cross product of the angular velocity vector and the unit vectors /, J, and K. The pro- 
gramming steps and the displayed output are 

(Cl) ENTERMATRIX( 3, 1 ) ; 


ROW 1 COLUMN 1 D[1 ,1 ] ; 
ROW 2 COLUMN 1 D[I,2]; 
ROW 3 COLUMN 1 D[I , 3] ; 
MATRIX-ENTERED 

(Dl) 


[ 

[ D 

[ 


1,2 


[ D 

[ 


1,3 


] 

] 

] 

] 

] 

] 


(C2) ENTERMATRIX( 3, 3) ; 


ROW 

1 

COLUMN 

1 

0; 

ROW 

1 

COLUMN 

2 

- P C 3 ] ; 

ROW 

1 

COLUMN 

3 

P[2]; 

ROW 

2 

COLUMN 

1 

P[3] ; 

ROW 

2 

COLUMN 

2 

0; 

ROW 

2 

COLUMN 

3 

-PCI]; 

ROW 

3 

COLUMN 

1 

-P[2]; 

ROW 

3 

COLUMN 

2 

PCI]; 

ROW 

3 

COLUMN 

3 

0; 
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MATRIX- ENTERED 
(D2) 


The product of these two matrices is 
(C3) (D2).(D1); 

(D3) 


[ 0 

[ 

E p. 


[ -p, 



£ ] 

■ p ’ ] 
o ] 
] 


■ P 2 D I ,3 “ P 3 D I ,2 =j 
= P 3 D I,1 " P 1 °I,3 jj 


The individual terms of this column vector can be evaluated for / = 1,2,3 by executing the fol- 
lowing program statement: 


(C4) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
EV(C[I,J]:R0W((D3),J))$ 


The evaluated terms can be printed out by using the now familiar display statement 


(C5) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO DISPLAY(C[I ,J])$ 


c 

= 

[D 

P - 

■ D 


P ] 

1 . 

,1 

[ 1: 

,3 2 

1 

,2 

3] 

c 

= 

[D 

P - 

- P 

D 

] 

1 ; 

,2 

[ 1, 

,1 3 

1 

1 

,3] 

c 

= 

[P 

D 

- D 


P ] 

1 

,3 

[ 1 

1,2 

1 

,1 

2] 

c 

— 

[P 

D 

- D 


P ] 

2 

,1 

[ 2 

2,3 

2 

,2 

3] 

C 

= 

[D 

P - 

■ P 

D 

] 

2 , 

,2 

[ 2 . 

,1 3 

1 

2 

: ,3] 


c 

: [P 

D 

- P 

D ] 

2,3 

[ 1 

2,2 

2 

2,1] 

C 

= [P 

D 

- P 

D ] 

3,1 

[ 2 

3,3 

3 

3,2] 

C 

= EP 

D 

- P 

D ] 

3,2 

[ 3 

3,1 

1 

3,3] 

C 

= [P 

D 

- P 

D ] 

3,3 

[ 1 

3,2 

2 

3,1] 
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The dependence of the direction cosines on the indices / and J and the time T can be shown by 
using the DEPENDENCIES statement. The use of this statement facilitates the formulation of the 
differential coefficients 


(C6) DEPENDENCIES (D( I ,J ,T) )$ 

It only remains to request that the differential coefficients of the direction cosines DCjj with 
respect to the time Tbe added to the coefficients C/j and displayed as follows: 


(C7 ) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO 
DC[I,J]:C[I,J]+DIFF(D[I,J],T)$ 

(C8) FOR 1:1 THRU 3 DO FOR J:1 THRU 3 DO D I SPLAY ( DC [ I ,J])$ 


DC 

= [D 


P ■ 

- D 


P ] 

+ 

d_ 

D 

= 0 

U 

[ 1 

,3 

2 

1, 

2 

3] 


dT 

1,1 


DC 

= [D 


P ■ 

- P 

D 

] 

+ 

d_ 

D 

= 0 

1,2 

[ 1 

,1 

3 

1 

1 

,3] 


dT 

1,2 


DC 

= [P 

D 


■ D 


P ] 

+ 

d_ 

D 

= 0 

1,3 

[ 1 

1 

,2 

1, 

1 

2] 


dT 

1,3 


DC 

= 'P 

D 


- D 


P ] 

+ 

d_ 

D 

= 0 

2,1 

; 2 

2 

,3 

2, 

2 

3] 


dT 

2,1 


DC 

= [D 


P - 

- P 

D 

] 

+ 

d_ 

D 

= 0 

2,2 

[ 2 . 

,1 

3 

1 

2 

,3] 


dT 

2,2 


DC 

= [P 

D 

_ 

• P 

D 

] 

+ 

d_ 

D 

= 0 

2,3 

[ 1 

2 

,2 

2 

2 

,1] 


dT 

2,3 


DC 

= 'P 

D 

_ 

■ P 

D 

] 

+ 

d 

D 

= 0 

3,1 

; 2 

3 

,3 

3 

3 

,2] 


dT 

3,1 


DC 

= [p 

D 

_ 

- P 

D 

1 

+ 

d_ 

D 

= 0 

3,2 

[ 3 

3 

,1 

1 

3 

,3! 


dT 

3,2 


DC 

= [p 

D 

_ 

• P 

D 

] 

+ 

d_ 

D 

= 0 

3,3 

[ 1 

3 

,2 

2 

3 

,1] 


dT 

3,3 
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This concludes the formulation of the simplified aeronautical model considered. The formula- 
tion gave rise to 18 equations: 3 force equations; 3 moment equations; 9 direction cosine equations 
to determine the spatial orientation of the vehicle; and 3 equations to determine the geographical 
location of the vehicle relative to an Earth-fixed reference frame. It is seen that the technique of 
symbolic mathematical computation, as implemented by the MACSYMA system, can be used to 
facilitate the formulation of complex mathematical models of physical systems and reduce the 
errors to which human operators are prone. The versatility and simplicity of the system make it 
attractive to programmers and nonprogrammers alike. Moreover, as already noted, the capability of 
working interactively enhances the utility of the system by permitting the user to modify the for- 
mulation as he proceeds. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, California 94035, January 12, 1979 
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